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The  lattice  Boltzmann  equation  (LBE)  method  is  a kinetics-based  approach  for  fluid 
flow  computations,  and  it  is  amenable  to  parallel  computing.  Compared  to  the  well- 
established  Navier-Stokes  (NS)  approaches,  critical  issues  remain  with  the  LBE  method, 
noticeably  flexible  spatial  resolution,  boundary  treatments,  and  dispersion  and  relaxation 
time  mode.  Those  issues  are  addressed  in  this  dissertation  with  improved  practice 
presented.  At  the  formulation  level,  both  the  single-relaxation-time  (SRT)  and  multiple- 
relaxation-time  (MRT)  models  are  analyzed.  The  SRT  model  involves  no  artificial 
parameters,  with  a constant  relaxation  time  regulating  the  physical  value  of  fluid 
viscosity.  The  MRT  model  allows  different  relaxation  time  scales  for  different  variables. 
Computational  assessment  shows  that  the  MRT  model  has  advantages  over  the  SRT 
model  in  maintaining  stability,  reducing  the  oscillation,  and  improving  the  convergence 
rate  in  the  computation. 
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A multi-block  method  is  developed  for  both  the  SRT  and  MRT  model  to  facilitate 
flexible  spatial  resolutions  according  to  the  flow  structures.  The  formulae  for 
information  exchange  at  the  interface  between  coarse  and  fine  grids  are  derived  to  ensure 
the  mass  and  momentum  conservation  while  maintaining  the  second-order  accuracy.  A 
customized  time  matching  between  coarse  and  fine  grids  is  also  presented  to  ensure 
smooth  exchange  information.  Results  show  that  the  multi-block  method  can  greatly 
increase  the  computational  efficiency  of  the  LBE  method  without  losing  the  accuracy. 

Two  methods  of  force  evaluation  in  LBE  are  examined:  one  based  on  stress 
integration  on  the  solid  boundary  and  the  other  momentum  exchange  between  fluid  and 
solid.  The  momentum  exchange  method  is  found  to  be  simpler  to  implement  while  the 
integration  of  stress  requires  evaluation  of  the  detailed  surface  geometry  and 
extrapolation  of  stress-related  variables  to  the  same  surface.  The  momentum  exchange 
method  performs  better  overall. 

Improved  treatments  for  both  opened  and  solid  boundaries  are  presented.  It  is 
demonstrated  that  these  treatments  are  simpler  than  those  existing  ones  while  offering 
second  order  accuracy.  Furthermore,  the  issues  associated  with  moving  solid  boundaries 
are  addressed.  The  proposed  technique  can  substantially  reduce  the  specious  pressure 
fluctuation  as  the  boundary  crosses  the  grid  lines. 
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CHAPTER  1 
INTRODUCTION 

Recently,  there  has  been  much  progress  in  developing  the  method  of  the  lattice 
Boltzmann  equation  (LBE)  (McNamara  and  Zanetti  1988,  Higuera  and  Jimenez  1989, 
Chen  et  al.  1992,  Qian  et  al.  1992)  as  an  alternative,  computational  technique  for  solving 
complex  fluid  dynamic  systems  (Benzi  et  al.  1992,  Chen  and  Doolen  1998).  Adopting 
the  macroscopic  method  for  computational  fluid  dynamics  (CFD),  the  macroscopic 
variables  of  interest,  such  as  velocity  u and  pressure  p,  are  usually  obtained  by  solving 
the  Navier-Stokes  (NS)  equations  (e.g.,  Peyret  and  Taylor  1983,  Fletcher  1988).  In  the 
LBE  approach,  one  solves  the  kinetic  equation  for  the  particle  velocity  distribution 
function  f{x,  £t),  where  £ is  the  particle  velocity  vector,  x is  the  spatial  position  vector, 
and  t is  the  time.  The  macroscopic  quantities  (such  as  mass  density  p and  momentum 
density  pu)  can  then  be  obtained  by  evaluating  the  hydrodynamic  moments  of  the 
distribution  function/  This  approach  was  first  proposed  by  Frisch  et  al.  (1986),  with  the 
additional  theoretical  foundation  established  in  the  subsequent  papers,  notably 
McNamara  and  Zanetti  (1988),  Higuera  and  Jimenez  (1989),  Koelman  (1991),  and  Qian 
etal.  (1992). 

A popular  kinetic  model  adopted  in  the  literature  is  the  single  relaxation  time 
approximation,  the  so-called  Bhatnagar-Gross-Krook  (BGK)  model  (Bhatnagar  et  al. 
1954): 

/+  £-v  /■ -!(/-/«)  (1.1) 

ot  A 


1 
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where  /cqy  is  the  equilibrium  distribution  function  (the  Maxwell-Boltzmann  distribution 
function),  and  A is  the  relaxation  time. 

To  solve  for  / numerically,  Eq.  (1.1)  is  first  discretized  in  the  velocity  space  using  a 
finite  set  of  velocities  {%a}  without  affecting  the  conservation  laws  (Chen  and  Doolen 
1998,  He  and  Luo  1997a  , 1997b), 

%-+4-v/„  = ~C4-^-’)  (1.2) 

ot  A 

In  the  above  equation,/^*,  t ) =j{x,  £,a,  t)  is  the  distribution  function  associated  with 
the  a-th  discrete  velocity  and  is  the  corresponding  equilibrium  distribution 

function.  The  9-velocity  square  lattice  model,  which  is  often  referred  to  as  the  D2Q9 
model  (Figure  1.1)  has  been  successfully  used  for  simulating  2-D  flows.  In  the  D2Q9 
model,  we  use  ea  to  denote  the  discrete  velocity  set  and  we  have 

eo  = 0, 

ea=c(cos((a  - l)zr / 4),  sin((<ar  - l)zr / 4))  for  a=l,  3,  5,  7, 

ea=j2c(cos((a-\)x/4),  sin((a- 1)^/4))  fora=2,4,6,8  (1.3) 

where  c = 5x/S /,  8x  and  5 1 are  the  lattice  constant  and  the  time  step  size,  respectively. 
The  equilibrium  distribution  for  D2Q9  model  is  of  the  form 

faq)  = Pwa[  1+-^-  eau+  -?—(ea-u)2  - uu ] (1.4) 

c2  2c4  2c 2 

where  wa  is  the  weighting  factor  given  by 


4/9, 

a = 0 

1/9, 

a - 1,3, 5, 7 

[1/36, 

a = 2, 4, 6, 8. 

With  the  discretized  velocity  space,  the  density  and  momentum  fluxes  can  be  evaluated 


as 
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P = Z/«=Z/“  (1.6) 

k=0  k= 0 

and 

P"  = Z ' eafa  = Z ( 1 -7) 

*=1  *=1 

The  speed  of  sound  in  this  model  is  cs  =c/^3  and  the  equation  of  state  is  that  of  an  ideal 
gas, 

p = pc]  (1.8) 

Equation  (1.2)  can  be  further  discretized  in  space  and  time.  The  completely 

discretized  form  of  Eq.  (1.1) , with  the  time  step  5t  and  space  step  eaSt,  is: 

fa(xt  + eadt,  t + &) -fa(x„  t)=--[fa(Xi,t)-f^\x„t)\  (1.9) 

T 

where  r =A/8t,  and  x,  is  a point  in  the  discretized  physical  space.  The  above  equation  is 
the  discrete  lattice  Boltzmann  equation  (McNamara  and  Zanetti  1988,  Higuera  and 
Jimenez  1989,  Chen  et  al.  1992)  with  BGK  approximation  (Bhatnagar  et  al.  1954).  Eq. 
(1.9)  is  usually  solved  in  the  following  two  steps: 


collision  step: 

7«  <*„  t + ay =/„  (at,,  o -h/«  (v  o - fP'  (*,. »)] 

T 

(1.1  Oa) 

streaming  step: 

fdxj  + eaa,  t + a)  =fa  (xi,  t + st) 

(1.10b) 

where  ~ represents  the  post-collision  state. 

The  viscosity  in  the  NS  equation  derived  from  Eq.  (1.8)  is 

v = (r -\/ 2)c]5t  (1.11) 

This  result  makes  formally  the  LBGK  scheme  a second  order  method  for  solving 
incompressible  flows  (He  and  Luo  1997a,  1997b).  The  positivity  of  the  viscosity 
requires  that  r>  Vi. 
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It  is  noted  that  the  collision  step  is  completely  local  and  the  streaming  step  takes  very 
little  computational  effort.  Eq.  (1.9)  is  explicit,  easy  to  implement,  and  straightforward  to 
parallelize. 

These  inherent  advantages  of  the  LBE  method  can  be  best  maintained  with  the  use  of  a 
regular  lattice  structure  (such  as  a square  lattice  or  hexagonal  lattice)  with  uniform 
spacing.  A challenge  of  using  the  uniform  grid  is  the  need  for  higher  resolution  in  high 
gradient  region  such  as  near  wall  high  Reynolds  number  flows,  and  where  reduced 
resolution  is  acceptable,  such  as  the  outer  boundary  far  away  from  the  body  for  external 
flows.  In  order  to  use  the  regularly  spaced  lattice  while  developing  the  capability  to  place 
the  outer  boundary  far  away,  it  is  desirable  to  divide  the  computational  domain  into  a 
number  of  grid  blocks  so  that  within  each  block  uniform  lattice  spacing  can  be  used. 
Such  a multi-block  approach  has  been  actively  employed  in  the  Navier-Stokes  (NS) 
equation  methods  with  both  Cartesian  and  curvilinear  coordinates  (Rai,  1985,  Steger, 
1991,  Shyy  et  el.,  1994).  A similar  multi-block  approach  in  LBE  will  be  desirable  also. 

Like  in  any  other  fluid  flow  computations,  the  numerical  boundary  condition  is  a very 
important  issue  in  the  LBE  method.  Two  types  of  boundaries  are  considered  for,  namely, 
the  solid  wall  and  the  open  boundary.  Unlike  solving  the  NS  equations  where  the 
macroscopic  variables,  their  derivatives,  or  a well  established  constraint  (such  as  the  mass 
continuity  for  pressure  distributions)  can  often  be  explicitly  specified  at  boundary  (Shyy 
et  al.,  1994),  in  the  LBE  method,  these  conditions  need  to  be  converted  into  distribution 
function/  Due  to  the  use  of  the  Cartesian  grids,  the  boundary  of  a practically  interesting 
geometry  will  often  intersect  with  the  grid  lines  irregularly.  To  obtain  an  accurate 
numerical  result,  it  is  important  to  resolve  the  geometric  details  in  such  situations.  The 
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bounce-back  scheme  is  a standard  boundary  condition.  It  is  intuitively  derived  from 
lattice  gas  automata  and  has  been  extensively  applied  in  lattice  Boltzmann  equation 
simulations.  In  this  scheme,  the  particles  hit  the  boundary  and  bounces  back  opposite  to 
its  incoming  direction  as  illustrated  in  Figure  1.2. 

Finally,  the  presence  of  moving  boundaries  is  of  substantial  interest  in  many  important 
practical  problems  (Shyy  et  el.,  1996).  The  treatment  of  the  moving  boundary  can  add 
substantial  complexities  and  burden  to  a given  method.  In  the  context  of  the  LBE 
method,  most  moving  boundaries  work  has  been  reported  to  address  the  particle 
suspension  in  fluids.  An  overall  review  for  the  application  of  lattice  Boltzmann  method 
to  simulation  of  particle-fluid  suspension  has  been  given  by  Ladd  and  Verberg  (2001). 

The  present  dissertation  is  structured  as  follows.  Chapter  2 will  review  in  detail  the 
development  of  the  theoretical  foundation  of  the  LBE  method  and  the  numerical 
implication.  Chapter  3 will  present  the  concepts  related  to  the  multi-block  techniques, 
with  special  attention  paid  to  ensure  the  conservation  of  mass  and  momentum  fluxes 
across  the  interface  between  blocks.  Next,  in  Chapter  4,  methods  for  evaluating  force 
(drag  and  lift)  on  a body  are  discussed.  In  Chapter  5,  a multi-relaxation-time  LBE  model 
is  introduced,  and  it  is  shown  computationally  that  this  model  can  reduce  the  spatial 
oscillations  observed  in  the  numerical  high  Reynolds  number  flow  solution.  In  Chapter  6, 
suitable  boundary  treatment  strategies  will  be  discussed  for  both  irregularly  shaped  solid 
walls  and  the  moving  boundary  problem  as  well  as  for  open  boundaries.  Finally, 
conclusions  and  future  work  are  discussed  in  Chapter  7. 
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Figure  1.1  A 2-D,  9-velocity  lattice  (D2Q9)  model. 


fluid 
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Figure  1.2  Bounce-back  inlet  condition  gives  f~(t)  = fa(t),  where  e-=-ea, 
and  ~ denotes  the  post-collision  state. 


CHAPTER  2 

BACKGROUND  OF  THE  LATTICE  BOLTZMANN  EQUATION  (LBE)  METHOD 

2.1  Lattice  Gas  Automata 

Broadwell  (1964)  presented  a kinetic  equation  with  single  particle  speed  to  simulate 
the  flow  involving  a shock  wave.  In  his  model,  time  and  space  are  continuous.  Hardy  et 
al.  (1976)  proposed  a fully  discretized  kinetic  model  on  square  lattice.  However,  due  to 
the  lack  of  isotropy  of  the  lattice,  this  model  failed  to  recover  the  NS  equation.  Frish, 
Hasslacher,  Pomeau  (1986)  presented  a 2-D  kinetic  model  based  on  the  hexagonal  lattice. 
In  this  model,  which  is  often  referred  to  as  the  FHP  model,  the  isotropy  condition  is 
enforced  and  for  the  first  time  the  NS  equation  is  successfully  recovered.  In  the  FHP 
model,  the  Boolean  variables  (0  or  1 ) are  used  to  describe  the  particle  occupation. 

The  construct  of  the  kinetic  model  involves  two  steps.  First,  a correct  symmetrical 
lattice  must  be  designed.  Second,  suitable  rules  of  evolution  must  be  established.  There 
are  two  alternatives:  lattice  gas  automata  (LGA)  and  lattice  Boltzmann  equation 
(McNamara  and  Zanetti  1988).  In  LGA,  Boolean  variables  are  used.  At  each  node,  a set 
of  Boolean  variables  na(x,t)(cc  = 1 ,...,M)  is  used  to  describe  the  particle  state,  where  M 
is  the  number  of  particle  velocities.  The  number  n represents  the  state  of  particle.  When 
na=  1,  there  is  a particle  in  the  ath-direction  and  when  na=  0 there  is  no  particle  in  the 

ath-direction.  The  evolution  equation  of  LGA  is  as  follows: 

na  (*  + ea  J + 1)  = na  (*»0  + («(*>0)  (2-1) 

where  ea  is  the  particle  velocities  and  Qa  is  the  collision  operator.  Evolution  involves 

the  following  two  steps: 
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Streaming:  particle  moves  to  the  nearest  neighboring  nodes  along  its  velocity 
direction. 

Collision:  when  arriving  at  nodes,  particles  collide  with  each  other  and  change  their 
velocities  and  directions  according  to  collision  rules. 

An  example  of  the  2-D  FHP  kinetic  model  is  shown  in  Figure  2.1  to  illustrate  the 
computational  procedures.  In  the  2-D  FHP  model,  the  hexagonal  lattice  is  used.  At  each 
node,  there  are  six  links  to  its  nearest  neighboring  nodes.  Starting  from  an  initial  state, 
the  particles  move  to  the  neighboring  nodes  with  unit  velocity  in  unit  time.  After  arriving 
at  the  new  location,  the  particles  collide  with  each  other  and  change  their  directions. 
Figure  2.2  shows  two  possible  outcomes  of  the  head-on  collision  of  two  particles  in  the 
FHP  model. 

The  construction  of  the  collision  rule  is  a crucial  part  in  the  LGA.  In  the  collision 
step,  the  mass,  momentum,  and  energy  conservation  need  to  be  maintained  while  no 
spurious  invariance  should  be  preserved  when  the  conservation  laws  are  enforced.  For 
example,  the  head-on  collision  can  results  in  two  possible  scenarios  as  shown  in  Figure 
2.2,  but  because  of  the  lack  of  redistribution  particles  along  any  part  of  the  opposite 
direction,  these  collision  rules  alone  will  generate  spurious  invariance.  The  problem  can 
be  solved  by  using  triple  collision  rules  or  putting  a rest  particle  at  every  node  (Firsch  et 
al.  1986). 

The  advantages  of  LGA  are  obvious.  The  evolution  rules  are  very  simple.  Only 
Boolean  operation  is  needed  in  the  computation.  Some  special  machine  can  be  designed 
and  built  to  take  the  advantage  of  these  characteristics.  The  drawback  is  also  obvious. 
The  macroscopic  variables  need  to  be  obtained  by  averaging  over  a large  number  of 
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possible  combinations  of  random  initial  condition  of  particle  states  over  a long  time 
(Benzi  et  al.  1992).  The  process  will  inevitably  produce  statistical  noise.  When  high 
order  physical  quantities  such  as  vorticity  or  stress,  which  involves  the  derivatives  of 
macroscopic  variables,  are  needed,  a much  lower  level  of  noise  is  needed. 

2.2  Lattice  Boltzmann  Equation 

2.2.1  Development  of  LGA 

Historically  the  LBE  method  was  derived  from  LGA.  In  LBE,  the  Boolean  variables 
in  the  LGA  are  replaced  by  real  variables,  fa  =<na  >,  in  which  < > represents  the 

ensemble  average.  Her e fa  is  the  particle  distribution  function.  Because  of  the  use  of  real 
variables,  the  collision  rules  need  to  be  redefined  in  LBE.  The  discrete  kinetic  equation 
for  particle  distribution  in  LBE  is 

fa  (x  + eaAt>t  + AO  = fa  (x,/)  + Qa  (/(x,0)  (2.2) 

where  Qa(/(x,/))is  the  collision  operator  which  represents  the  rate  of  change  of  fa  due 
to  collision.  At  is  the  time  step,  and  ea  is  the  discretized  velocity  in  phase  space  and 
|e„|  = A^/ . Eq.  (2.2)  requires  that  space  be  discretized  in  a way  that  is  consistent  with 

the  kinetic  equation.  The  coordinates  of  the  nearest  neighbor  points  around  x are  x+ea.  It 
is  noted  that  Qa  depends  on  the  local  distribution  function  and  the  lattice  structure. 

The  macroscopic  variables  such  as  density,  p,  and  momentum  density,  pu,  are 
calculated  from  the  moments  of  particle  distribution  function,/^: 

m m 

P = P‘*  = YJfaea  (2-3) 

a= 0 a=l 

where  m+1  is  the  number  of  discretized  velocity  in  the  phase  space.  The  collision 
operator  must  satisfy  the  conservation  of  mass  and  momentum  at  each  node. 
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m m 


I«„=o  Xn„e„=o 


(2.4) 


Since  At  physically  is  the  time  required  for  a particle  to  move  from  one  node  to  another 
and  is  much  smaller  than  the  time  scale  that  characterizes  the  variation  of  physical 
variables,  and  At  can  be  taken  as  a small  parameter  of  order  in  £ in  a perturbation 
analysis.  Using  a Taylor  series  expansion  in  space  and  time  in  Eq.  (2.2),  the  following 
continuous  form  of  the  kinetic  equation,  which  is  accurate  up  to  second  order  in  £,  can  be 
obtained. 


Assuming  that  the  distribution  function  f is  close  to  its  value  at  the  equilibrium  state, 
fa  can  be  expanded  as 


where  /a(I)is  the  leading  order  perturbation  that  accounts  for  the  major  part  of  the  non- 
equilibrium part  of fa.  Substituting  Eq.  (2.6)  into  the  collision  operator  Qa  in  Eq.  (2.5) 
and  performing  the  Taylor  series  expansion,  one  obtains 


(2.5) 


(2.6) 


(2.7) 


Substituting  Eqs.  (2. 6-2. 7)  into  (2.5),  we  find. 


na{fM)  = 0 


(2.8) 


This  leads  to  the  linearized  collision  operator: 


(2.9) 
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where  M aj} 


_ an,  (/<-’) 

at, 


the  coefficient  of  the  collision  matrix  (Higuera  and  Jimenez 


1989).  The  following  definitions  for  Map  are  given  by  Benzi  et  al.  (1992  ): 

(i)  The  [ Map\  is  a symmetrical  matrix. 

(ii)  The  general  element  Map  only  depends  on  the  angle  between  the  directions  of  ea  and 


ep. 

(iii)  Collision  conserves  mass  and  momentum,  so  that 

m m 

SM^=0’  0=1,... ,m  (2.10) 

a=0  a = 0 


The  matrix  elements,  Map,  are  numerical  parameters,  which  can  be  changed  at  will.  If  we 
further  assume  that  the  local  particle  distribution  relaxes  to  an  equilibrium  state  with  the 
same  relaxation  time  r,  it  can  be  obtained  that 

Map  =-~Sap  (2.11) 

T 

Subsequently  the  lattice  Boltzmann  equation  with  BGK  collision  (hereinafter  referred  to 
as  LBGK  equation)  is  obtained: 

/.  (*  + A»,  t + At)  = /„  (x,r)  - - (J,  - /“ ) (2.12) 

r 

The  BGK  collision  term  was  first  proposed  by  Bhatnagar,  Gross,  and  Krook  (1954). 

2.2.2  Derivation  of  LBE  from  the  Boltzmann  Equation 

Without  the  external  force,  the  Boltzmann  equation  with  BGK  collision  is  as  follows 
(Bhatnagar  et  al.  1954): 

%■+  f-v/  =-!(/-/<■»>) 

at  A 


(2.13) 
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where  / = f(x,^,t)  is  the  particle  distribution  function  in  continuum  phase  space,  £is 


the  particle  velocity,  and  feq  is  the  Maxwell-Boltzmann  equilibrium  distribution 
function: 


/ 


eq 


P 

(2  ttRT)DI2 


exp(- 


( Z-u)2 
2 RT  ' 


(2.14) 


where  T is  the  temperature,  R is  ordinary  gas  constant,  D is  the  spatial  dimension,  and  u 
is  the  flow  velocity.  The  macroscopic  variables  are  obtained  by  taking  various  moments 
of  the  distribution  function: 

p = [fdc  pu  = (2.15) 

Equation  (2.13)  is  a partial  differential  equation  and  we  can  use  any  known  numerical 
methods  to  solve  it.  After  discretizing  Eq.  (2.13)  in  phase  space  and  replacing  £with  e , if 
we  apply  the  first  order  finite  difference  scheme  in  time,  the  first  order  upwind 
discretization  for  the  convective  term  gives 

fa(x,t  + At)  = fa  (x,  0 - d[  fa  (x  + A x,t  + At)  - fa  (x,  t + AO] 


(2.16) 

where  d = Atea  / Ax,  and  At  and  Ax  are  the  time  step  and  grid  size,  respectively. 
Choosing  d=  1 and  denoting  r = XI  At , Eq.  (2.16)  becomes  standard  LBE, 

+ + = (2.17) 

T 


Although  only  the  first  order  schemes  are  used  to  obtain  the  discretized  differential 
equation,  as  shown  by  Sterling  and  Chen  (1996),  the  discretization  error  has  a special 
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form  which  can  be  included  into  the  viscous  term,  and  the  solution  to  Eq.  (2.17)  gives 
second  order  accuracy  for  the  macroscopic  variables  in  time  and  space. 

2.3  From  the  LBE  to  the  Navier-Stokes  (NS)  Equations 
In  order  to  recover  the  NS  equations  from  the  LBE  scheme,  we  begin  with  the  most 
commonly  used  single-relaxation-time  (SRT),  LBGK  model,  expressed  in  Eq.  (2.17). 
Performing  the  Taylor  series  expansion  in  time  and  space,  we  obtain 


At 


%!_ 

dt 


+ Ate 


ak 


dfa  , (A/)' 
dx.  2 


d2fa 


dt 


2 


52fa 

dtdxL 


+ eakecm 


d2fa 

dxkdxn 


+ -(/«  = 0 (2.18) 
r 

In  order  to  derive  the  NS  equations  from  LBE,  the  Chapman-Enskog  expansion  is  used. 
In  essence,  it  is  a standard  multi-scale  expansion  (Kervokian  and  Cole,  1980),  with  time 
and  space  being  expanded  as 

t = £t\  -t-  £~tj  + ' " X = £T|  + X j + • • • 


— = £ + £~ 

dt  dt. 


dt , 


d d 

= £ 

dx  dx. 


(2.19) 


and  the  particle  distribution  function  fa  expanded  similarly  as, 

fa=fJ0)+£ti')+£2f?)+O(£3) 

Substituting  equations  (2.19-2.20)  into  Eq.  (2.18)  results  in 

(0)  (0) 


i(/,T -/„“)+* 

T 


df  w df  w 1 

At^-  + Ateak^—  + -fc 


0) 


dt. 


dx 


\k 


(2.20) 


At 
V - 


3/®  3/a(1> 


dt. 


dt-. 


" ak 


dx 


\k 


+ 


(A  tr 


2 W0) 


d / dL 

lSL v-Jp  Ja 

2 ^ ^ ak 


dt 


dt,  dx 


+ eakecm 


\k 


svr 

3*1*  *ln 


y 
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s2~C=  0 

T 


Arranging  terms  in  ascending  orders  in  e,  we  obtain 

0(8°):  .C=faeq) 


0(e  ): 


/a0)  = -*A/ 


dfa 


(0) 


06 


■ + e 


(0) 


ak 


dx 


\k 


0(s2): 


/j2)=  ~rAt 


SfP  , . dfa 

+ eojt 


(1) 


0/,  0/, 


0r 


u 


-r 


(At)2 


2 r( 0) 


2 WO) 


0/ 


2 


ay,1 

0t,0X„ 


+ eaitean 


a2/j0) 

a***,* 


Using  Eq.  (2.23),  Eq.  (2.24)  can  be  simplified  to 


0f(O)  1 

f^  = -rAt^-+^-r)At 

otj  2 


dt , 0r 


li- 


lt is  noted  that  the  equilibrium  distribution  should  satisfy  following  constraints: 


P = Hfa 


(eq) 


o= I 


m 

Pu  = YJfaCq)ea 


a=l 


From  equations  (2.22)  and  (2.25)  we  find 


f.  =/„w+^',+£!/iu+o(£J) 


Equations  (2.25-2.26)  lead  to  following  results: 


Z/«(1)  =Z/J2)  = o 

o= 0 o=0 


Z famea  =Z  famea  =0 


o=0 


o=0 


(2.21) 

(2.22) 

(2.23) 

(2.24) 

(2.25) 

(2.26) 

(2.27) 

(2.28) 


(2.29) 
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Summing  Eq.  (2.23)  over  the  phase  space,  the  continuity  equation  is  recovered  to  the  first 
order  in  a. 


dp_+dpu± 
dt,  dx,. 


= 0 


(2.30) 


Multiplying  Eq.  (2.23)  with  the  discrete  particle  velocity  e ^ and  summing  over  «,  the 


momentum  equations  to  the  first  order  in  a are  recovered  as 


dt 


z 


+ > e^e 


Va 


(0) 


an  ak 


1 /= 1 


dx 


= 0 


(2.31) 


\k 


In  order  to  obtain  the  NS  equations  we  need  to  use  the  expression  for  f to  the  second 
order  of  a Summing  over  Eq.  (2.25)  and  enforcing  the  conditions  given  by  Eqs.  (2.28- 
2.29)  lead  to 

dp 


dt , 


= 0 


(2.32) 


Multiplying  Eq.  (2.25)  with  and  summing  over  a , we  obtain 
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Combining  equations  (2.30)  and  (2.32),  to  the  second  order  of  a,  the  continuity  equation 
equations  is  recovered: 


dP  , dpUp  _ Q 

dt  dxp 


(2.34) 


Combining  equations  (2.31)  and  (2.33),  to  the  second  order  of  a,  the  momentum 
equation  becomes 
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Changing  the  notation  of  k in  Eq.  (2.23)  to  m and  substituting  it  into  (2.35),  we  obtain  the 
momentum  equations  by  integrating  over  the  equilibrium  distribution  of f. 
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Using  the  definition  of  fa]  = \ckcJ(0)dVc  = pukun  + P5km  , one  obtains 


a=l 


dpun  d(pukun)  dp  ,1 


dt 


dx 


u 


dx 


+ (--r)A/Z 


1/7 


a=l 


2 ^ (0) 


£ , £ 
ak  cm 


dlfa 


dt\dxxk 


+ e ,e  e 

ak  an  am 


dxudx \m  J 


(2.37) 


In  Eq.  (2.35),  we  can  see  that  the  viscous  terms  depend  on  the  leading  order  term  of 
the  non-equilibrium  part  of  distribution  function,  while  equilibrium  part  only  determines 
momentum  flux  and  pressure  terms.  The  key  to  the  recovering  of  the  NS  equation  is  to 
find  an  equilibrium  distribution  function  in  the  discretized  phase  space. 

2.4  Two-Dimension  9-Velocity  Model  (D2Q9) 

The  9-velocity  (or  9-bit)  LBE  model  on  the  2-D  square  lattice  (Figure  1.1),  denoted  as 
the  D2Q9  model,  has  been  widely  used  for  simulating  2-D  flows.  For  athermal  fluids,  the 
equilibrium  distributions  are  defined  in  Eq.  (1.4).  To  recover  the  NS  equations  the 
following  identities  are  noted: 


8 c4 


'YjW<*ecmeaneakeaj  ~~^mnSkj  + S mk5 nj  + 5 mj5 „k) 


a=i 


(2.38) 


a d ^ 

2^Waean,ean^ake(,pUJ  =~  ( 


4 a2 


dX\ m dXk6  ff=l 


9 


(pun)  + 2 


dx^faxj 


puj  (2.39) 


Substituting  Eq.  (1.4)  into  the  second  order  derivative  terms  of  /(0)  in  Eq.  (2.36)  and 


using  the  identity  (2.39),  we  obtain 
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Substituting  Eq.  (2.40)  back  into  Eq.  (2.36),  we  obtain  the  momentum  equations: 
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where  the  dynamic  shear  and  bulk  viscosities  are 


c2.  1,. 

Ms = y(r--)Af 


(2.42) 


It  is  noted  that  r>l/2  is  required  to  obtain  a positive  viscosity.  In  the  incompressible 
flow  limit  we  have 
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dt  a dxa  p dxp 

The  differences  between  LBE  and  incompressible  NS  equation  are  as  follows: 

1.  LBE  is  a kinetic  equation;  NS  equations  are  macroscopic  equations. 

2.  The  convection  terms  in  LBE  are  linear;  NS  equations  have  nonlinear  convection 
terms. 

3.  LBE  is  a discretized  kinetic  equation;  NS  equations  can  take  either  integral  or 
differential  forms. 

4.  LBE  depends  on  the  grid  structure;  NS  is  in  vector  form  that  is  independent  on  the 
coordinate  and  grids. 

5.  Iterative  procedure  is  usually  employed  to  obtain  converged  solution  when  numerically 
solving  the  NS  equations.  LBE  does  not  employ  iterative  procedure. 
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6.  Generally  LBE  is  designed  for  time-dependent  simulations.  NS  equations  can  take 
either  steady  state  or  unsteady  form. 

2.5  Practical  Issues  in  the  Implementation  of  LBE  Method 
2.5.1  Boundary  Conditions 

Two  classes  of  boundaries  are  often  encountered  in  computational  fluid  dynamics: 
open  boundaries  and  the  solid  wall.  The  open  boundaries  include  lines  (or  planes)  of 
symmetry,  periodic  cross-sections,  infinity,  and  inlet  and  outlet.  At  these  boundaries, 
velocity  or  pressure  is  usually  specified  in  the  macroscopic  description  of  fluid  flow. 

A difficulty  of  the  LBE  method  is  that  the  boundary  conditions  for  the  distribution 
functions  / are  not  known.  One  must  construct  a suitable  condition  for^s  based  on  the 
macroscopic  flow  variables. 

At  the  symmetric  and  periodic  boundaries,  the  conditions  for/xS  can  be  given  exactly. 
At  the  outlet  of  a computational  domain,  the  condition  can  usually  be  given  by  simple 
extrapolation.  The  boundary  condition  for  the  velocity  at  a solid  wall  can  only  be  given 
approximately.  A popular  approach  is  to  employ  the  bounce-back  scheme  (Ladd  1994a, 
1994b,  Behrend  1995).  The  same  idea  for  the  solid  wall  has  been  directly  extended  at 
inlet  to  provide  the  required  inlet  condition  (Cornubert  et  al.  1991,  Ziegler  1993,  Behrend 
1995).  While  this  kind  of  treatment  is  an  easy  extension  of  the  solid  wall  condition,  it 
also  creates  a mechanism  for  the  pressure  waves  from  the  interior  of  the  computational 
domain  to  strike  on  the  inlet  boundary  and  reflect  back  into  the  computational  domain. 
Because  the  pressure  wave  propagation  is  typically  of  an  inviscid  nature,  depending  on 
the  nonlinearity  and  complexity  of  the  specific  problem,  this  may  provoke  computational 
instability  or  prevent  the  spatial  wiggles  from  being  dissipated  quickly  during  the  course 
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of  computation.  It  is  desirable  to  specify  a suitable  boundary  condition  for  / s that 
prevents  such  a feedback  of  the  pressure  wave  into  the  computational  domain. 

For  the  evolution  of  the  particle  distribution  function  / near  a solid  wall,  it  is  noted  that 
the  values  of / after  the  collision  at  the  nodes  inside  the  solid  wall  that  is  adjacent  to  the 
fluid  nodes  are  needed  in  order  to  complete  the  streaming  step.  Different  schemes  have 
been  proposed  to  obtain  post-collision  values  of / at  these  solid  nodes.  The  bounce-back 
method  is  very  simple,  efficient  and  easy  to  implement.  On  the  other  hand,  the  bounce- 
back  scheme  is  in  general  only  a first  order  treatment  (Ginzbourg  and  Alder  1994)  and 
can  reduce  the  overall  accuracy  of  the  solution.  Filippova  and  Hanel  (1998)  developed  a 
curved  boundary  condition  treatment  using  Taylor  series  expansion  in  both  space  and 
time  for  fa  near  the  wall.  This  boundary  condition  gives  a second  order  accurate 
treatment  for  a curved  solid  wall.  In  addition,  two  other  boundary  treatments  for  curved 
wall,  one  by  Mei  et  al.  (1999)  and  the  other  by  Bouzidi  et  al.  (2001)  have  also  been 
proposed.  All  of  those  methods  need  to  treat  the  boundary  condition  separately  for  A<V i 
and  zl>  Vi  in  which  zl  is  the  fraction  of  an  intersected  node-to-node  link  in  the  fluid 
region.  This  will  cause  non-uniformity  in  / when  A changes  from  bellow  Vi  to  above  Vi 
(or  vice  versa ) in  dealing  with  curved  boundaries.  A unified  scheme  for  curved  wall  is 
preferred. 

When  a solid  boundary  is  moving,  the  issue  is  more  complicated.  A special  concern 
emerges  when  a solid  boundary  moves  across  a node  so  that  it  leaves  behind  a new  fluid 
node  on  one  side  of  the  body  and  it  may  transform  a neighboring  fluid  node  into  a solid 
node  on  the  other  side  of  the  body.  The  specification  of / on  the  newly  generated  fluid 
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node  may  create  spurious  pressure  oscillations  which  can  degrade  the  solution  at  higher 
Reynolds  number. 

2.5.2  Reconciliation  of  Resolution  Requirement 

High  resolution  is  often  needed  in  the  high  gradient  region,  such  as  the  near  wall  flow 
of  a high  Reynolds  number  flow.  This  requires  a small  grid  size  near  the  wall.  For 
external  flows,  the  far  field  boundary  often  needs  to  be  placed  far  away  from  the  solid 
wall  while  the  resolution  can  be  lower.  This  immediately  creates  a difficulty  for  the 
standard  LBE  method,  since  it  inherently  requires  the  grid  size  to  be  constant.  Since  the 
grid  size  is  first  chosen  to  honor  the  high  gradient  region  near  the  wall,  it  requires 
enormously  large  numbers  of  grid  point  to  be  able  to  reach  the  far  field.  On  the  other 
hand,  since  the  gradients  of  physical  variables  are  small,  the  use  of  fine  grids  near  the  far 
field  is  a waste  of  computational  effort.  In  the  NS  equation-based  solvers,  this  issue  is 
easily  resolved  by  using  grid  stretching  techniques,  often  in  a curvilinear  coordinate 
system.  To  increase  the  numerical  efficiency  while  maintaining  accuracy  in  LBE,  non- 
uniform  grid  has  been  introduced.  He  et  al.  (1996)  developed  a time-dependent 
interpolation  scheme,  which  increases  grid  resolution  in  the  high  shear  rate  region.  Mei 
and  Shyy  (1999)  used  a curvilinear  coordinate  system  in  LBE  using  a finite  difference 
formulation. 

To  maintain  the  inherent  advantage  of  LBE,  one  prefers  to  employ  the  uniform  lattice. 
Thus  to  reconcile  the  different  requirements  of  grid  resolution,  it  is  desirable  to  divide 
computational  domain  into  different  regions.  In  each  region,  a different  lattice  space  can 
be  used.  In  practical  application,  there  are  two  different  treatments.  One  is  that  the  whole 
computational  domain  is  covered  with  coarse  grids.  The  patches  with  fine  grids  are 
placed  at  the  desired  regions.  The  other  is  that  there  is  only  one  layer  of  grids  at  each 
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region.  Such  a multi-block  approach  has  been  actively  employed  in  the  NS  equation- 
based  solution  methods  with  both  Cartesian  and  curvilinear  coordinates.  A similar  multi- 
block strategy  in  LBE  will  enhance  the  capability  of  the  LBE  method. 

2.5.3  Force  Evaluation 

The  evaluation  of  the  force  on  a body  in  a flow  field  is  a common  task  in  fluid 
mechanics  research  efforts.  This  is  also  the  case  for  computations  using  the  LBE  method. 
Different  force  evaluation  schemes,  including  momentum  exchange  (Behrend,  1995, 
Ladd,  1994a,  1994b)  and  integration  of  stress  (He  et  al.  1996),  have  been  used  to  evaluate 
the  fluid  dynamic  force  on  a curved  body  in  the  context  of  the  LBE  method.  The  method 
of  integration  of  stress  requires  knowledge  of  the  details  of  the  surface  geometry  and  can 
be  laborious  for  3-D  problems.  On  the  other  hand,  the  accuracy  of  the  momentum 
exchange  method  needs  to  be  examined  at  higher  Reynolds  numbers. 

2.5.4  Computational  Stability  and  Dispersion 

Computational  instability  and  dispersion  are  issues  of  concern  in  many  practical 
computation.  The  simplest  lattice  Boltzmann  equation  is  based  on  the  BGK  model,  which 
involves  the  single-relaxation-time  (SRT)  approximation.  In  contrast  with  the  SRT 
approximation,  the  multi-relaxation-time  (MRT)  lattice  Boltzmann  equation  method  has 
been  examined  by  d’Humieres  (1992)  and  more  recently  by  Lallemand  and  Luo  (2000). 
The  basic  idea  is  that  in  the  SRT  model,  the  bulk  and  shear  viscosity  are  both  determined 
by  the  same  relaxation  time.  The  MRT  model  attempts  to  relax  different  modes  with 
different  relaxation  time.  In  the  MRT  method  bulk  and  shear  viscosity  are  determined  by 
different  relaxation  time,  so  bulk  viscosity  can  be  adjusted  independently.  It  is  inevitable 
to  generate  pressure  wave  (acoustic  wave)  in  computation  because  of  initial  condition  or 
singularity.  Large  bulk  viscosity  is  helpful  to  damp  such  kind  of  wave.  Numerically,  this 
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could  reduce  substantially  the  amount  of  energy  associated  with  the  pressure  wave  and 
thus  improve  the  accuracy  of  the  unsteady  solution  or  the  convergence  toward  steady 
state. 

Lallemand  and  Luo  (2000)  performed  detailed  linear  analyses  on  the  dispersion, 
dissipation,  and  stability  characteristics  of  a generalized  lattice  Boltzmann  equation 
model  (d’Humieres,  1992).  It  was  found,  through  the  investigation  of  simple  flow 
problems,  that  the  multiple  relaxation  times  in  the  generalized  lattice  Boltzmann 
equations  leads  to  better  computational  stability  due  to  the  separation  of  the  relaxations  of 
the  various  kinetic  modes.  However,  their  theoretical  study  indicates  that  there  is  little 
difference  in  the  performances  on  the  high  wavenumber  range.  Many  fluid  flow 
problems  possess  singularities.  Inadequate  resolutions  in  such  regions  have  often 
resulted  in  spatial  oscillations  in  the  solution.  In  nonlinear  problems,  such  dispersion 
errors  tend  to  promote  numerical  instability.  A key  issue  is  whether  the  MRT  model  can 
improve  the  numerical  stability  via  reduce  the  dispersion  errors  (or  wiggles)  without 
increasing  the  computational  cost. 

2.6  Scope  of  the  Research 

A multi-block  method  will  be  developed  to  reconcile  the  vastly  different  requirement 
on  the  resolution  in  different  flow  regions.  Derivations  for  the  interface  information 
exchange  will  be  presented  to  highlight  issues  associated  with  both  formal  order  of 
accuracy  and  physical  conservation  laws.  The  performance  of  the  present  multi-block 
method  will  be  assessed  by  examining  the  solution  accuracy  in  various  flow  problems 
and  by  examining  the  conservation  of  mass,  momentum,  and  stresses  across  the  block 


interface. 
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The  force  evaluation  schemes  using  momentum  exchange  for  curved  boundary  will 
be  examined  in  flow  problems  involving  simple  and  complex  geometry.  The  momentum 
exchange  method  will  be  shown  to  offer  better  accuracy  and  to  be  much  more  robust  to 
use  than  the  stress  integration  method. 

The  detailed  comparison  of  computational  performance  by  using  the  MRT  and  SRT 
models  will  be  presented  for  several  flow  problems  possessing  either  flow  singularity  or 
geometric  singularity.  A multi-block  method  in  conjunction  with  MRT  model  will  also 
be  presented  in  which  the  interface  information  exchange  is  specifically  derived  for  the 
MRT  model. 

Finally,  various  issues  related  to  the  computational  boundary  conditions  are 
examined.  The  interaction  between  the  inlet  boundary  and  the  fluid  nodes  in  the  interior 
of  the  computational  domain  will  be  studied  along  with  a strategy  to  reduce  such 
interactions.  For  solid  wall  boundary  condition,  a unified  expression  for  treating  the  solid 
boundary  will  be  developed  and  assessed.  Special  attention  will  also  be  paid  to  further 
improve  the  boundary  condition  for  moving  walls. 


Figure  2.1  A hexagonal  lattice 
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Figure  2.2  Two  possible  results  of  a head-on  collision  in  lattice  gas 


CHAPTER  3 

DEVELOPMENT  OF  MULTI-BLOCK  METHODS  IN  THE  LBE  METHOD 

3.1  Introduction 

Accuracy,  stability,  and  efficiency  determine  whether  the  application  of  a particular 
scheme  is  successful  or  not.  In  the  LBE  method,  the  limitation  to  the  numerical 
efficiency  is  that  the  discretization  of  LBE  was  constrained  on  a special  uniform  lattice. 
A challenge  of  the  uniform  grid  is  to  offer  high  resolution  near  a solid  body  and  to  place 
the  outer  boundary  far  away  from  the  body  without  wasting  the  grid  resolution  elsewhere. 
In  order  to  use  the  regularly  spaced  lattice  while  developing  the  capability  to  place  the 
outer  boundary  far  away,  it  is  desirable  to  divide  the  computational  domain  into  a number 
of  grid  blocks  so  that  within  each  block  uniform  lattice  spacing  can  be  used.  Again,  such 
a multi-block  approach  has  been  actively  employed  in  the  Navier-Stokes  equation 
methods  with  both  Cartesian  and  curvilinear  coordinates. 

To  enhance  the  numerical  efficiency  and  accuracy,  the  non-uniform  grid  has  been 
introduced  into  the  LBE.  He  et  al.  (1996)  developed  a time-dependent  interpolation 
scheme,  which  increases  grid  resolution  in  the  high  shear  rate  region.  Mei  and  Shyy 
(1999)  used  the  curvilinear  coordinate  system  in  LBE  using  finite  difference  formulation. 
Bouzidi  et  al.  (2001)  presented  LBE  on  the  two-dimensional  rectangular  grid.  All  BGK 
scheme  based  methods  need  some  kinds  of  interpolations  during  the  computation  which 
may  increase  numerical  viscosity  and  reduce  accuracy. 

Filippova  and  Hanel  (1998)  proposed  a local  refinement  method  where  the  post- 
collision distribution  function  passing  the  interface  was  rescaled  to  satisfy  the 
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conservation  of  mass  and  momentum  and  continuity  of  stresses  across  the  interface 
joining  two  grid  systems.  Patches  of  fine  grids  are  used  in  certain  regions,  for  examples, 
around  a solid  body.  Values  of  the  distribution  functions  on  the  coarse  grid,  which  are 
coming  from  regions  of  finer  patches,  with  a large  gradient  of  hydrodynamic  variables, 
are  calculated  with  second-order  interpolations  in  space  and  time  in  the  boundary  nodes 
of  the  fine  grid. 

In  the  presently  developed  method  the  flow  field  is  divided  into  blocks.  In  each  block, 
the  grid  spacing  is  uniform  with  desired  resolution.  An  accurate  interface  treatment 
between  neighboring  blocks  is  derived  to  ensure  the  continuity  of  mass,  momentum,  and 
stresses  across  the  interface.  Furthermore,  a systematic  effort  has  been  made  to  evaluate 
the  performance  of  the  present  multi-block  method,  including  both  accuracy  and 
efficiency  aspects.  Several  test  cases  have  been  employed.  A lid-driven  cavity  flow  is 
computed  using  a single  block  with  uniform  grid  and  the  present  multi-block  method. 
The  results  are  compared  with  published  benchmark  results.  A channel  flow  with  a 
parabolic  velocity  profile  at  the  inlet  over  an  asymmetrically  placed  cylinder  at  /?e=T00 
(based  on  the  average  incoming  velocity)  is  computed  next  using  the  multi-block  method. 
Finally,  flow  over  NACA0012  airfoil  at  /?e=500-5000  is  computed.  The  present  study 
shows  that  the  multi-block  strategy  can  greatly  improve  the  computational  efficiency  of 
the  LBE  method. 

After  the  completion  of  the  present  work,  two  papers  were  published:  one  by  Lin  and 
Lai  (2000),  the  other  by  Kandhai  et  al.  (2000).  In  Lin  and  Lai’s  method,  the  coarse  base 
grid  covers  the  entire  physical  domain,  and  the  finer  grid  blocks  are  placed  at  regions 
where  local  grid  refinement  is  desirable.  The  simulation  is  first  carried  out  on  the  base 
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grid  level  at  a smaller  relaxation  time,  allowing  a rapid  propagation  of  boundary 
information  throughout  the  entire  domain.  At  a later  time,  the  computation  of  the  fine 
grid  variables  is  initiated.  The  dependent  variables  on  both  grid  levels  are,  then, 
advanced  in  time  simultaneously  with  the  fine  grid  boundary  conditions  obtained  from 
the  base  grid  solution  at  the  grid  interface.  The  method  of  Kandhai  et  al.  (2000)  is  based 
on  multiple  nested  lattices  with  increasing  resolution.  The  discrete  velocity  Boltzmann 
equation  is  solved  numerically  on  each  sub-lattice  and  interpolation  between  the 
interfaces  is  carried  out  in  order  to  couple  the  computations  in  different  blocks. 

Compared  with  the  above-cited  references,  in  the  present  method,  the  different  grid 
size  blocks  are  not  overlapped  between  each  other,  and  blocks  are  connected  only 
thorough  interface.  Also  details  related  to  the  conservative  properties  and  block-to-block 
coupling  are  investigated  more  directly. 

3.2  Basics  of  the  Multi-Block  Strategy  in  the  LBE  Method 
To  illustrate  the  basic  idea,  a two-block  system  (a  coarse  and  a fine,  as  shown  in 
Figure  3.1)  is  considered  in  the  derivation  for  the  interfacial  information  exchange.  The 
ratio  of  the  lattice  spacing  between  the  two-grid  system  is 

m = dxc/&Cf  (3.1) 

For  a given  lattice  size  dx,  the  viscosity  of  the  fluid  is 

v = (2t  -\)Sxc  / 6 (3.2) 

In  order  to  keep  a consistent  viscosity,  and  thus  Re,  in  the  entire  flow  field  involving 
different  lattice  sizes,  the  relation  between  relaxation  times,  x/ , on  the  fine  grid,  and  rc, 
on  the  coarse  grid,  must  obey  the  following  rule: 

1 , K 

rf=-  + m(Tc--) 


(3.3) 
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To  keep  the  variables  and  their  derivatives  continuous  across  an  interface  between  two 
different  grids,  a consistent  and  accurate  relationship  for  the  probability  density  function 
in  the  neighboring  grid  blocks  must  be  developed.  The  following  summarizes  the  key 
elements  in  the  derivation  for  the  information  exchange  across  the  interface. 

It  is  noted  that 


/„(*.  o =/"(*  o + fr>(*.  i) 


(may) 


(3.4) 


where  f(neq)  is  the  non-equilibrium  part  of  the  distribution  function  based  on  which  the 


deviatoric  stresses  are  evaluated.  The  collision  step  gives 


1 


fa(x„  t + St)=  (1  —)/„(*„  + t) 

T T 


(3.5) 


Substituting  Eq.  (3.4)  into  Eq.  (3.5)  leads  to 
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fa(X„  t + St ) - (1  t)  + fr\xit  0]  t) 

r L J r 


/“(*„  0+— o 

T 


(3.6) 


Denoting  the  coarse-grid  quantities  with  the  superscript  c and  fine-grid  quantities  with  the 
superscript/  the  post-collision  step  gives 


?(c)  _ f(eq,c)  , ^ r(neq,c) 

J a J a ' J a 

r 


(3.7) 


Similarly, 


~ T r 1 

f(f)  _ f(eqj)  . L f f(neqj) 

J a Ja  ' Ja 


(3.8) 


Since  the  velocity  and  density  must  be  continuous  across  the  interface  between  the  two 
grids,  from  Eq.(1.6),  it  is  seen  that 

f(eq,c)  _ r(eqj) 

J a J a 


(3.9) 
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To  maintain  continuity  in  the  deviatoric  stresses,  in  the  2-D  case, 


it  is  required  that 


i 1 

r = (1— L)Y  /■<"■»> (e  e --e  e 8) 

lJ  v ^ ' / jJ  a V ai  a j ^ a a ij  J 


Qf=l 


(1~)  /j”*'1’  - (1-T1-)  /j”*-71 

2r.  2r , 


(3.10) 


(3.11) 


or 


f(aneq'c)  =m^~  f{aneqJ) 

Tf 

Substituting  Eq.  (3.12)  into  Eq.  (3.7)  one  obtains 


= fimJ) 


Using  Eqs.  (3.8-3.10),  the  above  becomes 


?(c)  = I 

J a J a 


T - 1 T 


m- 


r -1 


(3.12) 


(3.13) 


Tf  Tf 


if)  _ /•(««./)  i 


r/~1 


In  transferring  the  data  from  the  coarse  grids  to  the  fine  grids,  one  similarly  obtains 


(/)  _ J f 


fin=  /•< 

-/a  7 a 


r,  -1 


w(rc  - 1) 


\flc)~fleq’c)\ 


(3.15) 


Equations  (3.14-3.15)  were  first  given  by  Filippova  and  Hanel  (1998).  On  the  interface 
between  two  blocks,  there  are  m values  of  /„(/)  needed  for  each  f(eq'c)  and  f(c).  Thus, 

spatial  and  temporal  interpolation  procedures  for  the  values  of  f(aeq'c)  and  /je)  on  the 
fine-grid  lattice  are  used  to  complete  the  evaluation  of  f(ar) . 

3.3  Interface  Structure  and  Computational  Procedure 


The  typical  interface  structure  is  shown  in  Figure  3.1.  The  line  MN  is  the  fine  block 


boundary,  while  the  line  AB  is  the  coarse  block  boundary.  The  coarse  block  boundary  is 
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in  the  interior  of  the  fine  block,  and  the  fine  block  boundary  is  in  the  interior  of  the  coarse 
block.  This  arrangement  of  the  interface  is  convenient  for  the  information  exchange 
between  two  neighboring  blocks.  For  example,  grid  Q is  an  interior  lattice  node  of  the 
coarse  block.  After  the  collision  step,  the  values  of  incoming  distribution  functions 

f2(tn+[,XD),  f3(tn+',XE),  and  f4(t"+',XF)  from  boundary  nodes  D,  E,  and  F, 

respectively,  are  needed  in  order  to  obtain  f2(tn+\X0),  f3(tn+l,Xa),  and  f4(tn+',XQ)  at 

the  end  of  steaming  step,  since  other  components  of  fa(tn+\ X0),(a  = 1,5, 6, 7, 8)  are 

obtained  from  advecting  the  neighboring  post-collision  values  of  fa  in  the  interior  nodes 

of  the  coarse  block.  For  the  same  reason,  the  fine  block  boundary  MN  is  located  in  the 
interior  of  the  coarse  block.  However,  on  the  fine  block  boundary  MN , there  is  no 
information  on  the  nodes  denoted  by  the  solid  symbol  • in  Figure  3.1;  it  must  be  obtained 
through  spatial  interpolation  based  on  the  information  at  the  nodes  denoted  by  the  open 
symbol  ° on  the  line  MN . To  eliminate  the  possibility  of  spatial  asymmetry  caused  by 
interpolations,  a symmetric,  cubic  spline  fitting  is  used  for  spatial  interpolation  of  fa  on 
the  fine  block  boundary, 

f(x)  = al+blx  + clx2+dlxi,  x(_,<x<x(.,  i = \,...,n  (3.16) 

where  the  constants  ( al,bi,ci,dj ) are  determined  by  using  the  continuity  of  the  nodal 

conditions  of  /,/ ,/  and  suitable  end  conditions  (such  as  zero  second  derivative  for/). 
We  found  that  it  is  very  important  to  maintain  the  spatial  symmetry  in  the  interpolation 
along  the  interface. 
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Figure  3.2  shows  the  flow  chart  of  the  computational  procedure  for  the  multi-block 
calculation.  It  is  noted  that  in  addition  to  the  spatial  interpolation,  there  is  a need  for 

temporal  interpolation  on  all  nodes  at  the  fine  block  boundary  MN  in  order  to  obtain 
fa(t  2 , MN) . Here  a 3-point  Lagrangian  formula  is  used: 

(3-17> 

*-i  7-r*i-*y 

j*k 

3.4  Results  and  Discussions 

Several  well-documented  flow  problems  are  selected  to  highlight  the  performance  of 
the  present  method.  In  all  cases,  the  boundary  condition  for  fa  in  the  solid  region  near  a 
wall  is  obtained  using  the  formulations  given  by  Mei  et  al.  (1999)  for  a curved  geometry. 

3.4.1  Lid-Driven  Cavity  Flow 

The  lid-drive  cavity  flow  has  been  extensively  used  as  a benchmark  solution  to  test  the 
accuracy  of  a numerical  method.  In  this  flow,  two  singular  points  at  the  upper  comers  of 
the  lid  require  high  resolution  to  obtain  satisfactory  stress  distribution  near  the  corner 
points.  To  assess  the  LBE  results,  the  benchmark  solutions  by  Ghia  et  al.  (1982)  are  used 
for  comparison. 

The  computations  are  carried  out  using  (i)  a single-block  with  uniform  lattice 
(129x129)  with  the  walls  placed  halfway  between  lattices,  and  (ii)  a multi-block  whose 
layout  is  shown  in  Figure  3.3.  For  the  multi-block  case,  in  the  two  upper  corner  regions, 
the  grid  resolution  is  increased  by  a factor  of  4.  For  Re=100,  the  relaxation  time  is  rc= 
0.56  for  the  coarse-grid  block  and  Tj  =0.74  for  the  fine-grid  block.  The  upper  wall 
velocity  is  17=0.0156.  The  initial  condition  for  density  is  unity  and  that  for  velocity  is 
zero.  The  streamlines  shown  in  Figure  3.4  are  obtained  from  the  single  block  solution 
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and  the  pattern  is  not  discernable  from  those  of  the  multi-block  solution.  The  positions  of 
the  centers  of  the  primary  vortices  are  (0.6154,  0.7391)  and  (0.6172,  0.7390)  for  uniform 
grid  and  multi-block  solutions  respectively,  compared  well  with  the  value  (0.6172, 
0.7344)  from  Ghia  et  al.  (1982).  The  u-  and  v-components  of  the  velocity  along  the 
vertical  line  and  horizontal  line  through  the  geometry  center  are  shown  in  Figure  3.5a  and 
Figure  3.5b,  respectively.  It  is  seen  that  while  the  single  block  method  with  129x129 
lattices  can  capture  most  of  physical  variables  satisfactorily,  the  multi-block  method  can 
improve  the  quality  of  the  solution.  Figure  3.6  shows  the  pressure  contour  from  the 
single-block  computation.  Because  of  the  singularity  at  the  upper  comers,  the  density 
contours  exhibit  noticeable  oscillations  due  to  the  insufficient  resolution  near 
singularities.  Figure  3.7  shows  the  pressure  contours  obtained  from  the  multi-block 
solution.  Significant  improvement  in  the  smoothness  of  the  solution  for  the  pressure  field 
over  that  of  the  single  block  solution  is  observed. 

In  an  NS  solver  for  incompressible  flows,  because  of  the  decoupling  of 
thermodynamic  pressure  and  velocity  field,  it  is  crucial  to  maintain  the  mass  conservation 
of  the  entire  flow  domain.  This  issue  becomes  more  critical  when  the  multi-block 
method  is  used  (Shyy  1997,  Shyy  et  al.  1999).  Also  for  incompressible  flows,  the 
pressure  is  arbitrary  up  to  a constant.  Hence  coupling  the  pressure  term  while 
maintaining  the  mass  flux  conservation  is  very  important.  Generally  speaking,  it  is 
difficult  to  maintain  simultaneously  the  continuity  of  mass,  momentum,  and  stresses 
across  the  interface  between  neighboring  blocks  because  interpolations  are  applied  to 
each  dependent  variable.  In  the  present  multi-block  LBE  method,  the  continuities  of 
mass  and  stresses  are  ensured  through  the  use  of  Eqs.  (3.14-3.15).  The  most  important 
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point  is  that  interpolations  are  only  applied  to  fi ’s  along  the  interface  and  this 

automatically  ensures  the  consistency  in  the  transfer  of  various  flux  terms  across  the 
interface. 

To  validate  the  above  arguments,  pressure,  shear  stress,  mass  flux  and  momentum  flux 
near  the  block  interfaces  are  examined  next.  Figure  3.8  shows  a local,  enlarged  view  of 
the  pressure  contour  around  an  interface  comer  point  indicated  by  the  circle  in  Figure  3.7. 
Clearly,  the  pressure  is  rather  smooth  across  the  interface  with  the  coarse-to-fine  grid  size 
ratio  of  m= 4.  Figure  3.9-Figure  3.1 1 show  the  contours  of  shear  stress,  mass  flux,  and 
momentum  flux  pux  . It  is  seen  that  these  physical  quantities  are  all  smooth  across  the 
interface. 

To  demonstrate  this  issue  more  clearly,  macroscopic  physical  quantities  on  one  part  of 
the  interface  (i.e.  line  A-B  in  Figure  3.3)  are  plotted  in  Figure  3.12-Figure  3.17.  After  the 
streaming  step  there  is  no  physical  value  on  the  interface  for  the  fine  grid.  Here  we  use  a 
second  order  extrapolation  to  obtain  the  fine  grid  value  on  the  interface.  Figure  3.12- 
Figure  3.15  show  that  mass  and  momentum  fluxes  match  very  well  between  the  fine-  and 
coarse-grid.  Figure  3.16  shows  the  shear  stress  profile.  In  most  part  of  the  interface  the 
fine-  and  coarse-  grid  solutions  agree  very  well  with  each  other.  The  discrepancy  appears 
near  the  upper  wall.  It  is  noted  that  for  in  the  fine-grid  blocks,  the  top  moving  wall  is 
located  half-way  between  two  horizontal,  fine-grid  lattices  with  a distance  of  A/Sxf 
=0.5  Sxf.  In  the  coarse-grid  block,  the  distance  between  the  wall  to  the  nearest  lattice  in 
the  fluid  region  is  Acdxc  =0.55c/  =0.5(Stc  /4=0.125  Sxc  for  m= 4.  This  mismatch  (A/*zlc) 
will  result  in  different  errors  in  the  boundary  condition  for  /J’s.  This  subsequently  affects 
the  accuracy  of  the  shear  stress  near  the  corner  of  the  block  interface  and  the  wall.  The 
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same  problem  also  appears  in  Figure  3.17  for  pressure.  The  behavior  of  the  solution, 
however,  can  be  easily  improved  by  using  a set  of  uniform-sized  fine  grids  for  the  fluid 
region  near  the  entire  upper  wall. 

3.4.2  Channel  Flow  over  an  Asymmetrically  Placed  Cylinder  at  Re=100 

Schafer  and  Turek  (1996)  reported  a study  of  a laminar  flow  over  a circular  cylinder 
placed  asymmetrically  inside  a channel.  The  cylinder  center  to  the  upper  wall  distance  is 

4.2  cylinder  radii,  and  the  cylinder  center  to  the  lower  wall  distance  is  4.0  radii.  In  the 
present  LBE  computation,  two  grid  sizes  are  used  with  the  coarse-to-fine  lattice  spacing 
ratio  m= 4 and  the  flow  domain  is  divided  into  five  blocks.  The  finer  grid  block  is  placed 
around  the  cylinder  as  shown  in  Figure  3.18.  The  radius  of  the  cylinder,  r,  equals  20 
lattice  spacing  in  the  fine  block.  The  total  number  of  lattices  in  coarse-grid  blocks  is 
8854,  and  the  fine  block  has  81  x 8 1=6561  lattices.  The  relaxation  times  for  the  coarse- 
and  fine-  grids  are  rc=  0.52  and  zy  =0.58,  respectively.  The  average  inlet  velocity  U is 
0.0666.  The  channel  inlet  has  a parabolic  velocity  and  is  located  at  4.0  radii  upstream  of 
the  cylinder  center.  A zeroth-order  extrapolation  for  fa  is  used  at  the  outlet.  The 
Reynolds  number  based  on  the  average  inlet  velocity  and  the  diameter  of  the  cylinder  is 
Re=  100. 

At  this  Reynolds  number,  the  flow  becomes  unsteady  and  periodic  vortex  shedding  is 
observed.  The  numerical  value  of  Strouhal  number  (St  = D/UT)  is  0.300  and  it  agrees 
very  well  with  the  range  of  values  (0.2995-0.305)  given  by  Schafer  and  Turek  (1996). 
Here  D is  the  diameter  of  cylinder  and  T is  the  peak-to-peak  period  of  the  lift  force  which 
is  500  in  lattice  unit  based  on  the  coarse  block.  An  instantaneous  streamline  plot  is 
shown  in  Figure  3.18  after  the  dynamically  periodic  solution  is  established.  The 
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variations  of  the  drag  and  lift  coefficients  are  shown  in  Figure  3.19.  It  is  noted  that 
CD{t ) has  two  peaks  (CDmax,  =3.23,  CDmax2  =3.22  ).  A closer  examination  of  C,(t) 

reveals  that  CL(t)  is  not  symmetric  with  respect  to  the  x-axis 
(CL max  = 1.01,  CLmm  =-1.03).  This  result  is  reasonable  because  the  flow  is  not 

symmetric  with  respect  to  the  horizontal  line  drawn  through  the  center  of  the  cylinder. 
The  mean  velocity  of  the  flow  passing  the  upper  region  of  cylinder  is  lower  than  that 
passing  the  lower  region  (Schafer  and  Turek,  1996).  Careful  examination  of  the 
computational  flow  filed  indicates  that  the  local  pressure  in  the  upper  region  is  higher 
than  that  in  the  lower  region  at  similar  stages  of  the  vortex  shedding.  There  is  no  report 
of  the  existence  of  two  peaks  of  Cn(t ) and  asymmetry  of  C,{t)  in  Schafer  and  Turek’s 

paper.  Only  the  ranges  of  CDmax  (3.22-3.24)  and  C,  max  (0.99-1.01)  were  given.  The 
present  results  for  CD  max  and  CL  max  are  well  within  those  ranges. 

3.4.3  Steady  Flow  over  the  NACA0012  Airfoil 

The  NACA  0012  airfoil  (Figure  3.20)  is  a popular  wing  model,  which  has  been  used 
extensively.  Flow  fields  at  Re=500,  1000,  2000,  and  5000  are  computed  with  the  multi- 
block LBE  scheme.  Figure  3.21  shows  the  entire  computational  domain  and  the 
schematic  of  the  multi-block  arrangement.  There  are  300  lattices  (grids)  along  the  chord 
in  the  finest  block.  The  largest  grid  size  ratio  between  neighboring  blocks  is  4.  At  the 
inlet,  upper,  and  lower  boundaries,  the  equilibrium  values  are  used  for  f s according  to 
Eq.  (1.4)  based  on  the  free-stream  velocity.  At  the  downstream  boundary  a zeroth  order 
extrapolation  for  f' s is  used. 
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Figure  3.22  shows  the  density  contour,  streamlines  and  velocity  vector  of  the 
converged  solution  at  Re= 2000  and  zero  angle  of  attack.  To  investigate  the  effect  of  grid 
resolution,  two  sets  of  grid  systems  are  used  for  the  flow  field  at  7?e=500:  a fine  grid 
system  and  a coarse  grid  system  (with  resolution  reduced  by  a factor  of  2 in  every  block). 
Figure  3.23  shows  the  velocity  profiles  at  {x-xle)!L= 0.06  where  L is  the  chord  length  and 
xie  is  the  location  of  the  leading  edge.  These  two  sets  of  velocity  profile  agree  well  with 
each  other,  although  the  fine  grid  solution  appears  to  have  smoother  u-component 
velocity  profile,  as  expected. 

Figure  3.24  compares  the  drag  coefficient  Q between  the  present  LBE  simulation  and 
those  calculated  from  Xfoil,  which  is  a coupled  inviscid  and  boundary  layer  flow  solver 
(Drela  and  Giles  1987).  It  can  be  seen  that  two  sets  of  results  agree  with  each  other  very 
well  for  the  range  of  Reynolds  numbers  investigated  in  this  study. 

It  is  also  noted  that  at  /?e=500,  the  present  value  of  Q=0.1761  compares  very  well 
with  the  results  reported  by  Lockard  et  al.  (2000):  Q=0.1762  obtained  using  a Navier- 
Stokes  equation-based  finite  difference  method,  and  C^=0.1717  using  Powerflow  code 
developed  by  EXA  Corporation,  which  is  based  on  the  lattice  Boltzmann  equation 
method.  In  addition,  the  present  simulation  for  the  symmetrical  flow  at  Re=500  gives  a 
lift  coefficient  of  |Q|  <6  x 10'14.  Lockard  reported  Cl  =1 . 1 5x1 0‘7  using  an  NS  equation- 
solver  and  Cl  =2.27x1  0"4  using  EXA’s  Powerflow  code.  This  suggests  that  the  present 
multi-block  code  preserve  the  symmetry  very  well. 

It  is  worth  pointing  out  that  there  is  a significant  saving  in  the  computational  cost 
using  the  multi-block  method  in  LBE  simulations.  There  are  three  different  sizes  of  grids 


used  for  the  NACA0012  airfoil  simulation.  There  are  1025x129=132225  fine  grids, 
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93200  intermediate  grids  with  m= 4,  and  139628  coarse  grids  (with  m= 8 in  reference  to 
the  finest  grids).  This  gives  a total  of  about  3.6xl0?  grids  in  the  entire  domain.  If  the  fine 
grid  system  is  used  in  the  entire  domain,  the  number  of  the  grids  would  be  NxxNy= 
5698x1153-  6.57xl06  which  is  18  times  more  than  in  the  multi-block  case.  This 
represents  a saving  of  18  times  in  the  memory.  Furthermore,  since  dt=5x=5y  in  the  LBE 
simulation,  one  time  step  in  the  coarsest  grid  system  (m= 8)  requires  2 time  steps  in  the 
intermediate  grid  blocks  and  8 steps  in  the  finest  grid  blocks.  The  ratio  of  the 
computational  efforts  required  to  carry  out  a single-block  simulations  to  that  for  a multi- 
block simulation  for  a given  period  of  physical  time  would  be 

6. 57xl06x8/(l 32225x8+93200x2+1 39628)  - 38. 

Clearly,  more  saving  can  be  achieved  if  more  blocks  of  different  sizes  are  used. 

3.5  Conclusions 

A multi-block  strategy  is  developed  for  the  lattice  Boltzmann  method.  The  interface 
condition  is  derived  to  ensure  the  mass  conservation  and  stress  continuity  between 
neighboring  blocks.  Favorable  computational  results  are  obtained  in  selected  test  cases. 
With  efficiency  aspects  greatly  improved,  there  is  a significant  potential  for  the  multi- 
block strategy  in  the  LBE  method  for  practical  flow  problems. 
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Figure  3.1  Interface  structure  between  two  blocks  of  different  lattice  spacing. 
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Figure  3.2  Flow  chart  of  the  computational  procedure  in  the  multi-block  method. 
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129*129  A 129*129 


Figure  3.3  Block  layout  for  a 2-D  cavity.  Lattice  spacing  is  reduced  by  a factor  of 
8 for  graphical  clarity. 


Figure  3.4  Streamlines  in  the  cavity  flow  at  Re=100. 
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Figure  3.5a  Comparison  of  u-velocity  along  the  vertical  line  through  geometric 
center. 


Figure  3.5b  Comparison  of  v-velocity  along  the  horizontal  line  through 
geometric  center. 

Figure  3.5  Comparison  of  velocity  between  present  results  and  those  by  Ghia  et  al. 
(1982). 
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Figure  3.6  Pressure  contours  in  the  cavity  flow  from  the  single-block  LBE 
simulation. 


Figure  3.7  Pressure  contours  in  the  cavity  from  multi-block  LBE  solution.  (For  the 
circled  region,  see  Figure  3.8) 
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Figure  3.8  Enlarged  view  of  pressure  contour  in  the  circled  region  in  Figure  3.7 
near  the  intersection  of  three  blocks.  The  figure  demonstrates  that  the  block 
interface  and  comer  are  well  handled. 


Figure  3.9  Shear  stress  contour.  Solid  and  dash  lines  represent  positive  and 
negtive  values,  respectrively. 
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Figure  3.10  Contour  of  x-component  mass  flux  pux.  Solid  and  dash  lines  represent 
positive  and  negtive  values,  respectrively. 


Figure  3.11  Contour  of  momentum  flux  in  the  x-direction  pu  ] . 
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Figure  3.12  The  x-component  of  the  mass  flux  pux  /( poll)  on  the  interface  AB.  In 
Figure  3.12-Figure  3.15,  po=  1 and  U=0.0156. 


Figure  3.13  The  x-component  of  the  mass  flux  puy  /(poU)  on  the  interface  AB. 
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Figure  3.14  The  x-component  of  the  momentum  flux  , pu2  / p0U2 , on  the  interface 
AB. 


Figure  3.15  The  y-component  of  the  momentum  flux,  puxuy  / p0U2 , on  the 
interface  AB. 
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Figure  3.16  Shear  stress  rv  /(/jU/H) on  the  interface  AB. 


Figure  3.17  Pressure  on  the  interface  AB. 


Figure  3.18  Instantaneous  streamlines  for  channel  flow  over  an  asymetrically 
placed  cylinder  at  Re=100. 
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Figure  3.19  Unsteady  drag  and  lift  coefficients  on  the  cylinder. 
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Figure  3.20  NACA  0012  airfoil. 


Figure  3.21  Block  and  lattice  layout  for  flow  over  NACA  0012.  The 
lattice  spacing  is  reduced  by  a factor  32  for  graphical  clarity. 
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Figure  3.22  Streamlines,  pressure  contour,  and  velocity  vectors  for  a uniform  flow 
over  NACA  0012  airfoil  at  Re=2000. 
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Figure  3.23  Grid-independence  test  of  the  velocity  profiles  near  the  leading  edge 
at  (x-xLe)/L=0.06  for  flow  over  NACA0012  airfoil  at  Re=500. 


Re 

Figure  3.24  Comparison  of  Cd  between  the  present  simulation  and  Xfoil 
calculation  as  a function  of  Re  for  flow  over  NACA0012  airfoil.  The  straight  line 
is  the  slope  according  to  the  laminar  boundary  layer  theory. 


CHAPTER  4 

FORCE  EVALUATION  ON  A SOLID  BODY  IN  THE  LBE  MTHOD 

4.1  Introduction 

4.1.1  Force  Evaluation  and  Related  Works 

In  spite  of  numerous  improvement  for  the  LBE  method  during  the  last  several  years, 
one  important  issue  that  has  not  been  systematically  studied  is  the  accurate  determination 
of  the  fluid  dynamic  force  involving  curved  boundaries.  Needless  to  say,  accurate 
evaluation  of  the  force  is  crucial  to  the  study  of  fluid  dynamics,  especially  in  fluid- 
structure  interaction.  Several  force  evaluation  schemes,  including  momentum  exchange 
(Behrend,  1995,  Ladd,  1994a)  and  integration  of  stress  (He  et  al.  1996),  have  been  used 
to  evaluate  the  fluid  dynamic  force  on  a curved  body  in  the  context  of  the  LBE  method. 

He  and  Doolen  ( 1 997)  evaluated  the  force  by  integrating  the  total  stress  on  the  surface 
of  the  cylinder  and  the  components  of  the  stress  tensor  were  obtained  by  taking  respective 
velocity  gradients.  Even  though  a body-fitted  grid  was  used,  an  extrapolation  was  needed 
to  obtain  the  stress  in  order  to  correct  the  half-grid  effect  due  to  the  bounce-back 
boundary  condition.  Filippova  and  Hanel  (1998)  developed  a second-order  accurate 
boundary  condition  for  curved  boundaries.  However,  the  fluid  dynamics  force  on  a 
circular  cylinder  asymmetrically  placed  in  a 2-D  channel  was  obtained  by  integrating  the 
pressure  and  deviatoric  stresses  on  the  surface  of  the  cylinder  by  extrapolating  from  the 
nearby  Cartesian  grids  to  the  solid  boundary  (Filippova  and  Hanel,  1998).  To  gain 
insight  into  the  method  of  surface  stress  integration,  it  is  instructive  to  examine  the 
variation  of  the  pressure  on  the  surface  of  a circular  cylinder  at  finite  Reynolds  number 
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obtained  using  LBE  method  for  flow  over  a column  of  cylinders  (Filippova  and  Hanel, 
1998).  Figure  4.1  shows  the  definition  of  boundary  points,  xw.  Figure  4.2  shows  the 
pressure  coefficient,  cp  = (p-p^Vf/i-vU2),  on  the  surface  obtained  by  using  2nd  order 
extrapolation  where  p*,  is  the  far  upstream  pressure.  Only  those  boundary  points,  xw, 
which  are  intersected  by  the  horizontal  or  vertical  lattice  velocity  vectors  ea,  (a=l,  3,  5, 
and  7)  are  included  in  Figure  4.2.  If  the  boundary  points  intersected  by  the  links  in  the 
directions  of  ej,  e4,  e$,  eg,  are  included,  the  variation  of  cp  would  be  more  noisy.  The 
components  of  the  deviatoric  stress  tensor  show  a similar  noisy  pattern.  It  is  not  clear 
how  the  noise  in  the  pressure  and  stresses  affect  the  accuracy  of  the  fluid  dynamic  force 
in  the  stress  integration  method.  While  the  programming  in  the  extrapolation  and 
integration  is  manageable  in  2-D  cases,  it  is  rather  tedious  in  3-D  cases. 

Instead  of  the  stress  integration  method,  Ladd  (1994)  used  the  momentum  exchange 
method  to  compute  the  fluid  force  on  a sphere  in  suspension  flow.  In  the  flow  simulation 
using  the  bounce-back  boundary  condition,  the  body  is  effectively  replaced  by  a series  of 
stairs.  Each  segment  on  the  surface  has  an  area  of  unity  for  a cubic  lattice.  The  force  on 
each  link  (half-way  between  two  lattices  at  xj  and  x*  =(xy  + eaSt ) in  which  xg  resides  in 
the  solid  region)  results  from  the  momentum  exchange  (per  unit  time)  between  two 
opposing  directions  of  the  neighboring  lattices, 

[eafa(xf’t)~eafa  (xb  = xf  + e a8t  ,t)]l  St  in  which  e~=-ea.  Whereas  the  momentum 

exchange  method  is  very  easy  to  implement  computationally,  its  applicability  and 
accuracy  for  a curved  boundary  have  not  been  systematically  studied. 

To  recapitulate,  there  are  two  major  problems  associated  with  the  method  of  surface 
stress  integration.  First  of  all,  the  components  of  stress  tensor  are  often  noisy  on  a curved 
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surface  due  to  limited  resolution  near  the  body  and  the  use  of  Cartesian  grids.  The 
accuracy  of  such  a method  has  not  been  addressed  in  the  literature.  Secondly,  the 
implementation  of  the  extrapolation  for  Cartesian  components  of  the  stress  tensor  to  the 
boundary  surface  and  the  integration  of  the  stresses  on  the  surface  of  a 3-D  geometry  are 
very  laborious  in  comparison  with  the  intrinsic  simplicity  of  the  lattice  Boltzmann 
simulations  for  flow  field. 

The  problems  associated  with  the  method  of  the  momentum  exchange  are  as  the 
following.  (1)  The  scheme  as  proposed  for  the  case  with  A=  V2  at  every  boundary 
intersection  xw.  Whether  this  scheme  can  be  applied  to  the  case  with  zfc  Vi,  when,  for 
example,  the  geometry  is  not  straight,  needs  to  be  investigated.  (2)  As  in  the  case  of 
stress  integration  method,  the  resolution  near  a solid  body  is  always  limited  and  the  near 
wall  flow  variables  can  be  noisy.  If  one  uses  the  momentum  exchange  method  to 
compute  the  total  force,  it  is  not  clear  what  is  the  adequate  resolution  to  obtain  reliable 
fluid  dynamic  force  on  a bluff  body  at  a given  (moderate)  value  of  Reynolds  number,  say, 
Re~0(  102). 

4.1.2  Scope  of  the  Present  Work 

In  what  follows,  two  methods  for  the  force  evaluation,  i.e. , the  stress  integration  and 
the  momentum  exchange  methods,  will  be  described  in  detail.  The  shear  and  normal 
stresses  on  the  wall  in  a pressure  driven  channel  flow  will  be  first  examined  to  assess  the 
suitability  of  the  momentum  exchange  method  when  A * Vi  and  analyze  the  errors 
incurred.  The  results  on  the  drag  force  for  flow  over  a column  of  circular  cylinders  using 
these  two  methods  will  be  subsequently  assessed  for  the  consistency.  The  drag 
coefficient  at  Re=100  will  be  compared  with  the  result  of  Fomberg  (1991)  obtained  by 


54 


using  a second  order  accurate  finite  difference  method  with  sufficient  grid  resolution. 
For  flow  over  an  asymmetrically  placed  cylinder  in  a channel  at  Re=100,  the  unsteady 
drag  and  lift  coefficients  match  well  with  literature  results.  The  unsteady  flow  over  The 
momentum  exchange  method  is  further  evaluated  for  3-D  fully  developed  pipe  flow  and 
for  a uniform  flow  over  a sheet  of  spheres  at  finite  Reynolds  number.  We  found  that  the 
simple  momentum  exchange  method  for  force  evaluation  gives  fairly  reliable  results  for 
2-D  and  3-D  flows. 

4.2  Methods  for  Force  Evaluation  in  the  LBE  Method 
4.2.1  Force  Evaluation  Based  on  Stress  Integration 

He  and  Doolen  (1997)  evaluated  the  force  by  integrating  the  total  stresses  on  the 
surface  of  the  cylinder, 

F = J n [-/?/  + pv(Vu  + (Vm)t)]  dA  (4.1) 

In  He  and  Doolen’ s study,  a body  fitted  coordinate  system  together  with  grid 
stretching  was  used  such  that  a large  number  of  grids  can  be  placed  near  the  body  to  yield 
reliable  velocity  gradient  V«.  In  general,  since  u is  not  the  primary  variable  in  the  LBE 
simulations  and  the  evaluation  of  u using  Eq.(1.7)  based  on  /„’ s suffers  the  loss  of 
accuracy  due  to  the  cancellation  of  two  close  numbers  in  f„ s,  the  evaluation  of  the 
derivative  Vh  will  result  in  further  degradation  of  the  accuracy.  Fillipova  et  al.  (1998) 
used  similar  integration  scheme  to  obtain  the  dynamic  force  on  the  body  for  the  force  on 
a circular  cylinder  except  that  the  deviatoric  stresses  were  evaluated  using  the  non- 
equilibrium part  of  the  particle  distribution  function  [see  Eq.  (4.3)  below].  However, 
since  Cartesian  grid  was  used,  the  stress  vectors  on  the  surface  of  the  body  (with  arbitrary 
A)  have  to  be  computed  through  an  extrapolation  procedure  based  upon  the  information 
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in  the  flow  field.  This  leads  to  further  loss  of  accuracy  for  finite  lattice  size  Sc  when  the 
shear  layer  near  the  wall  is  not  sufficiently  resolved. 

In  Eq.  (4.1),  the  pressure  p can  be  easily  evaluated  using  Eq.  (1.8).  For  D2Q9  and 

2 

D3Q19  models  cs  =1/3  so  that  p=pl3.  The  deviatoric  stress  for  2-D  incompressible  flow 

T,j  = pv  (V,  Uj  + V,  Ui ) (4.2) 

can  be  evaluated  using  the  non-equilibrium  part  of  the  distribution  function 

^ a= 1 ^ 

In  the  above,  f^neq)-  fa  - f^eq)  is  the  non-equilibrium  part  of  particle  distribution 
function /a. 

For  flow  over  a circular  cylinder,  a separate  set  of  surface  points  on  the  cylinder  can 
be  introduced  in  order  to  carry  out  the  numerical  integration  given  by  Eq.  (4.1).  The 
values  of  the  pressure  and  each  of  the  six  components  of  the  symmetric  deviatoric  stress 
tensor  on  the  surface  points  can  be  obtained  using  a second-order  extrapolation  scheme 
based  on  the  values  of  p and  zij  at  the  neighboring  fluid  lattices.  The  force  is  obtained  as 

F~  ^surface  ^ i~P  L £ (extrapolated  SA . (4.4) 


It  is  worth  commenting  here  that  for  the  2-D  flow  over  the  cylinder,  nearly  half  of  the 
entire  code  was  taken  up  by  the  above  force  evaluation  procedure.  Following  Fornberg 
(1991),  the  drag  coefficient  over  a circular  cylinder  of  radius  r is  defined  as 


CD  = 


1^1 

pU2r 


(4.5) 
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4.2.2  Force  Evaluation  Based  on  the  Momentum  Exchange 

In  order  to  employ  the  momentum  exchange  method  efficiently,  two  scalar  arrays,  w(i, 
j)  and  wh(i,  j)  are  introduced.  A value  of  0 is  assigned  to  w(i,  j)  for  the  lattice  site  (i,  j) 
that  are  occupied  by  fluid;  a value  of  1 is  assigned  to  w(i,j)  for  those  lattice  nodes  inside 
the  solid  body.  The  array  wt(i,  j)  is  set  to  zero  everywhere  except  for  those  boundary 
nodes,  Xb,  where  a value  of  1 is  assigned.  For  a given  lattice  direction  ea,  (a=  1,2,  ...,  8 in 
2-D),  its  opposite  direction  is  e~=  (- ea ) (see  Fig.  (4.1)).  For  a given  boundary  node  Xb 

inside  the  solid  region  with  w*(z,  j)=\  and  w(i,  /)=  1,  the  momentum  exchange  with  all 
possible  neighboring  fluid  nodes  over  a time  step  St=  1 is 

Nd 

I ea[fa(xb>t)  + fa(xb  +ea> t)][\~w(xb+ea)] 

a- 1 

where  Nj=S  for  2-D  and  Nj  =1 8 for  3-D  cases.  Simply  summing  the  contribution  over 
all  boundary  nodes  xh  belonging  to  the  body,  the  total  force  (acted  by  the  solid  body  on 
the  fluid)  is  obtained  as 

F=  X X ea  y a (xb . ' 0 + fa  (x„  + , 0]  [1  - w(xh  + ea  )]  (4.6) 

all  xh  a=\ 

In  the  momentum  exchange  method  the  force  F is  evaluated  after  the  collision  step  is 
carried  out  and  the  boundary  value  f~(xb,t ) has  been  evaluated  (Fillipova  et  al.,  1998, 
Mei  et  al,  1 999).  The  momentum  exchange  occurs  during  the  subsequent  streaming  step 
when  f-(xb,t)  and  fcfecf,  t)  move  to  xj  and  x^  respectively.  As  mentioned  in  the 

introductory  section,  the  effect  of  variable  A is  not  explicitly  included,  but  the  effect  of  A 
is  implicitly  taken  into  account  in  the  determination  of  f„(xb,t)  (Fillipova  et  al.,  1998, 
Mei  et  al,  1999).  The  applicability  of  Eq.  (4.6)  will  be  examined  and  validated. 
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Clearly,  the  force  is  proportional  to  the  number  of  boundary  nodes  x t,  in  the  above 
formulation  and  the  number  of  the  boundary  nodes  increase  linearly  with  the  size  of  the 
body  in  2-D  flow  in  the  LBE  method.  However,  since  the  force  is  normalized  by  pU2r  as 
in  Eq.  (4.5)  in  2-D,  the  drag  coefficient  Co  should  be  independent  of  r. 

4.3  Results  and  Discussions 

For  straight  walls,  there  is  no  doubt  that  Eq.  (4.1)  together  with  Eq.  (1.8)  for  pressure 
and  Eq.  (4.3)  for  zjj  gives  accurate  result  for  the  force  provided  that  /a’s  are  accurately 
obtained.  To  demonstrate  the  correctness  of  Eq.  (4.6)  based  on  the  momentum  exchange 
for  an  arbitrary  A,  consideration  is  given  first  to  the  pressure  driven  channel  flow  (Figure 
4.3)  for  which  exact  solutions  for  the  velocity  and  stresses  are  known.  The  values  of  the 
drag  on  the  column  of  circular  cylinder  at  100  for  Htr= 20  computed  using  these  two 
methods  are  then  compared  with  the  result  of  Fomberg  (1991).  The  dependence  of  the 
drag  on  the  radius  r in  the  momentum  exchange  method  is  examined  to  assess  the 
reliability  of  this  method. 

For  3-D  flows,  the  pressure  driven  flows  in  a circular  pipe  have  known  exact  solutions 
for  both  the  velocity  profile  and  wall  shear  stresses.  The  assessment  for  the  momentum 
exchange  method  in  3-D  flow  will  be  made  first  in  this  case.  Finally,  the  momentum 
exchange  method  will  be  evaluated  by  considering  the  drag  on  a sphere  due  to  a uniform 
flow  over  a sphere  in  a finite  domain.  The  details  for  flow  field  computation  can  be 
found  in  Mei  et  al.  (1999)  and  Mei  et  al.  (2000). 

4.3.1  Pressure-Driven  2-D  Channel  Flow 

In  the  case  of  the  channel  flow,  the  force  on  the  top  wall  (y=H)  at  a given  x location 
(i=Nx/ 2+1,  say)  can  be  evaluated  using  the  momentum  exchange  method  as  follows.  The 
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wall  is  located  between  j—Ny  and  Ny- 1 (Figure  4.3).  The  x-component  force  on  the  fluid 
at  the  wall  near  the  z'th  node  is 

F*  = t/6 (hj)+  /2 (/—  1 ’ 7— 1 )]  e6x  + [f%(j,j)+  74(z+1,;-1)]  e8x  (4.7a) 

Fy  = [?6 (f>j)  + /20'-l,7-l)]  e6y  + [7s  (Uj)+  TtO'+Ui'-l)]  e8y 

+[  fi  c uj ) + 73  O'.y-i )]  (4.7b) 

Since  <Sc=l,  Fx  and  Fy  are,  effectively,  the  total  shear  and  normal  stresses,  cr^  and  ayy, 
which  include  the  pressure  and  the  deviatoric  stresses,  on  the  fluid  element  at  y=H. 

Based  on  Eq.  (4.3),  the  deviatoric  component  of  the  fluid  shear  stresses  at  j=Ny- 1 (or 
y=Ny-3+A)  and  Ny- 2 (or  y=Ny-4+A)  can  be  exactly  evaluated  based  on  the  non- 
equilibrium part  of  the  distribution  functions  fa’s  in  the  flow  field  if  they  are  exactly 
given.  A linear  extrapolation  of  the  deviatoric  shear  stresses  to  y=H=  Ny-3+2A  yields 

ZZ  = Fy(j'=Ny-  l)  + A[  rxy(j=Ny- 1 ) - rxy(j=Ny-2J]  (4.8) 

where  the  superscript  “ neq ” denotes  that  the  value  is  based  on  using  fZq),  the 
subscript  w refers  to  the  value  on  the  wall.  The  deviatoric  normal  stress,  t"v^w  , can  be 

similarly  found.  In  a fully  developed  channel  flow,  the  normal  component  of  the 
deviatoric  stress  %,(y)  is  expected  to  be  zero  while  the  total  normal  stress  a}y(y)  is  equal 
to  the  pressure  (-/?).  It  needs  to  be  pointed  out  that  this  method  of  evaluating  fZ  given 

by  Eq.  (4.8)  for  2-D  channel  flow  is  equivalent  to  the  method  of  the  surface  stress 
integration  based  on  the  extrapolated  pressure  and  the  deviatoric  stresses  on  the  solid  wall 
except  that  no  numerical  integration  on  the  solid  surface  is  needed. 
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After  the  velocity  profile  ux(y)  is  obtained  from  fa  using  Eq.  (1.7),  the  shear  stress 
on  the  wall  can  also  be  calculated  using  the  near  wall  velocity  profile  as. 


M- 


du; 

dy 


| y=H  = pV 


,2  + A 0 - ux(j  = Ny-\) 
1 + A A 


-pv 


1 + A 


[ux(j=Ny-\)-ux(j=Ny-2)\  (4.9) 


In  the  above,  a linear  extrapolation  is  employed  to  evaluate  the  velocity  derivative 


— — \v=n  at  the  wall. 
dy  » 


Finally,  the  exact  solution  for  the  fluid  shear  stress  on  the  wall  (y=H)  is 

T™aw=~—H  with  H = Ny-  3+2A  (4.10) 

^’w  2 dx  y 

based  on  the  parabolic  velocity  profile  or  simple  control  volume  analysis.  This  exact 
result  can  be  used  to  assess  the  accuracy  of  the  aforementioned  method  for  the  force 
evaluation. 

In  the  LBE  simulations,  the  pressure  gradient  is  enforced  through  the  addition  of  an 
equivalent  body  force  (He  et  al.  1997,  Mei  et  al.  2000)  after  the  collision  step.  While  the 
velocity  field  given  by  the  LBE  solution  can  be  unique,  the  pressure  field  (thus  the 
density  field  p(x,  y))  can  only  be  unique  up  to  an  arbitrary  constant.  In  view  of  Eq.  (4.9), 
it  is  difficult  to  compare  the  stresses  for  different  cases  if  p(i,  j)  converges  to  different 
values  in  each  case.  To  circumvent  this  difficulty,  the  density  field  in  the  channel  flow 
simulation  is  normalized  by  p(i— 2,  j=Ny/2)  at  every  time  step.  This  normalization 
procedure  results  in  p(x,  y)=  1 throughout  the  entire  computational  domain.  It  is  also 
applied  to  the  3-D  flows  in  the  circular  pipe. 
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Table  4.1  compares  the  numerical  values  of  the  shear  stress  for  a typical  case  ( Ny=35 , 
dp/dx=  — 1 0"6,  and  r=  0.6)  based  on:  given  by  Eq.  (4.10),  Fx  given  by  Eq.  (4.7a), 

t7v\  given  by  Eq.  (4.8),  and  ^ \y=H  given  by  Eq.  (4.9).  Also  listed  is  the  comparison 

dy 

between  Fy  (Eq.  (4.7b))  and  -p.  All  computations  are  carried  out  with  double  precision 
accuracy. 

It  is  noted  that  is  identical  to  for  all  values  of  A.  Closer  examination  of 

the  shear  stress  profile  Tnxyq  (y)  using  Eq.  (4.3)  across  the  channel  reveals  that  zxyq(y)  is 

also  equal  to  the  exact  shear  stress  profile  z^act(y) , which  is  linear,  despite  the  errors  in 
the  velocity  profile  ux(y ) for  all  values  of  A.  A linear  extrapolation,  Eq.  (4.8),  for  a linear 
profile  therefore  gives  the  exact  wall  shear  stress.  Thus,  the  exactness  of  Txyqw  in  the  LBE 
simulation  of  channel  flow  indicates  the  reliability  of  the  LBE  solution  for  the  stress  field 
rnp  (x,  y)  using  Eq.  (4.3).  However,  as  Figure  4.2  indicated,  the  accuracy  of  integrating 

r"eq  (x,y)  to  obtain  the  fluid  dynamic  force  in  non-simple  geometries  is  not  clear;  this 
will  be  further  addressed  in  the  following  sections. 

For  0<  A <1,  the  normal  force  Fy  using  Eq.  (4.7b)  based  on  the  momentum  exchange 
method  agrees  exactly  with  the  pressure  on  the  wall.  This  is  a rather  special  quantity 
since  deviatoric  component  of  the  force  is  identically  zero.  Nevertheless,  the  method  of 
the  momentum  exchange  does  give  a reliable  value  for  the  normal  stress. 

For  the  shear  (tangential)  force,  it  is  observed  from  Table  1 that  for  fixed  dp/dx,  Fx 
does  not  change  as  A increases  from  0.01  to  0.99.  On  the  other  hand,  the  exact  result, 
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=——(Ny  - 3 + 2A) , increases  linearly  with  A.  Further  computations  were  carried 

y'  2 dx 

out  over  a range  of  Ny  (=35,  67,  99,  and  131)  and  r(=0.505,  0.51,  0.52,  0.6,  0.7,  0.8,  0.9, 
1.0,  1.2,  1.4,  and  1.6).  The  results  indicate  that  the  momentum  exchange  method  gives 
the  shear  stress  on  the  top  wall  as 

F*=\$r(Ny- 3 + |).  (4.11) 

2 ax  3 

That  is  Fx  is  not  sensitive  to  r and  A.  The  error  in  Fx  is  zero  when  A = 1/3.  The 

4 

absolute  error  attains  the  maximum  when  A = 1 which  gives  the  relative  error  of for 

3 H 

Fx.  Although  the  frequently  used  momentum  exchange  method  is  a natural  choice  for  the 
force  evaluation  in  conjunction  with  the  bounce-back  condition  for  A - 'A,  one  must  be 
aware  that  this  method  is  not  exact  and  the  error  in  the  force  evaluation  using  the 
momentum  exchange  method  depends  on  A and  resolution. 

Table  4.1  also  shows  that  for  the  shear  stress  based  on  taking  the  derivative  of  the 
velocity,  the  loss  of  accuracy  is  quite  significant  for  small  values  of  A (<0.05)  for  z=0.6. 
For  other  values  of  A (>0.3),  the  accuracy  is  comparable  with  that  of  Fx.  Flowever,  as 

shown  in  Figure  4.4,  the  accuracy  of  \y=H  based  on  the  near-wall  velocity 

dy 

derivative  deteriorates  as  the  relaxation  time  r increases  (from  0.51  to  1.6).  To  see  the 

du  i 

cause  of  the  increasing  error  in  fj. — — \y=H,  Figure  4.5  shows  dimensionless  wall  velocity, 

dy 

ujuc,  obtained  by  a 3-point  second-order  Lagrangian  extrapolation  of  the  near  wall  ux{y) 
as  a function  of  r.  The  increasing  slip,  uw,  on  the  wall  with  the  increasing  relaxation 
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time,  t,  was  also  observed  by  Noble  et  al.  (1995).  It  is  the  result  of  increasing  particle 
mean  free  path  that  causes  the  deviation  of  the  microscopic  based  solution  from  the 

dux  | 

continuum-based  solution.  It  is  clear  that  the  poor  performance  of  /u L =h  1S 

dy 

associated  with  the  increasing  error  in  the  near  wall  velocity  profile  as  r increases.  Since 
the  stress  tensor  r can  be  calculated  directly  from  fa  (Eq.  (4.3))  without  the  need  for 

taking  velocity  derivatives,  the  force  evaluation  method  based  on  the  evaluation  of  the 
velocity  gradient  in  the  form  of  Eq.  (4.2)  is  not  recommended. 

4.3.2  Steady  Uniform  Flow  over  a Column  of  Cylinders 

For  a uniform  flow  over  a column  of  circular  cylinders  of  radius  r with  center-to- 
center  distance  denoted  by  H (see  the  left  part  of  Figure  4.9  for  illustration),  symmetry 
conditions  for  fa’s  re  imposed  at  y=  ±H/2.  Most  of  the  details  of  flow  field  simulation 
can  be  found  in  Mei  et  al.  (1999).  The  Reynolds  Number  is  defined  by  the  diameter  of 
the  cylinder  d as  Re=Ud/v,  where  U is  the  uniform  velocity  at  the  inlet.  It  must  be  noted 
that  for  a consistent  determination  of  the  force,  the  upstream  boundary  must  be  placed  far 
upstream.  A shorter  distance  between  the  cylinder  and  the  boundary  will  result  in  higher 
drag.  In  this  study,  it  is  placed  at  about  20  radii  to  the  left  of  the  center  of  the  cylinder. 
Reducing  the  distance  between  boundary  and  the  cylinder  to  12.5  radii  while  keeping  the 
rest  of  the  computational  parameters  fixed  would  increase  the  drag  coefficient  by  about 
1.8%  at  Re=100.  The  downstream  boundary  is  located  about  25-30  radii  behind  the 
cylinder  to  allow  sufficient  wake  development.  The  simulation  is  terminated  when  the 
following  criterion  based  on  the  relative  L2-norm  error  in  the  fluid  region  Q is  satisfied. 


63 


{ Z [u(xhyj,t  + \)-u(Xj,yj,t)]2}112 

F - <, 

^2  ^ 1/7 

[ Z n2(Xj,yj,t  + 1)]1/2 

(yj,zk)eQ 


(4.12) 


In  this  case,  e=10'6  was  chosen  for  both  Re=10  and  100. 

Following  the  Fornberg  (1991),  the  drag  coefficient  over  a circular  cylinder  of  radius  r 
is  defined  as 


CD  = 


F ' 


pU2r 


(4.13) 


Figure  4.6a  compares  Co  obtained  from:  momentum  exchange  method,  surface  stress 
integration,  and  finite  difference  result  of  Fornberg  (1991)  using  a vorticity-stream 
function  formulation  at  Re=100,  Hlr=  20,  for  r ranging  from  2.8  to  13.2.  For  r>8,  both 
methods  of  momentum  exchange  and  the  stress  integration  give  satisfactory  results  for 
Cd  in  comparison  with  the  value  of  1.248  given  by  Fornberg  (1991).  The  small 
difference  in  Cd  could  be  due  to  the  fact  that  in  Fornberg’s  study,  the  computational 
domain  is  much  larger  in  downstream  direction-the  downstream  boundary  condition  is 
imposed  at  300  radii  behind  the  cylinder,  as  opposed  to  25-30  radii  in  present  study.  This 
adds  credence  to  the  validity  of  using  Eq.  (4.6)  for  evaluating  the  total  force  on  a body. 
The  values  of  Cd  from  the  momentum  exchange  method  have  a little  less  scattering  than 
that  from  the  stress  integration.  Accepting  an  error  of  less  than  5%,  the  reliable  data  for 
Cd  can  be  obtained,  using  the  momentum  exchange  method,  for  r> 5.  That  is,  10  lattices 
cross  the  diameter  of  the  cylinder  are  necessary  to  obtain  reliable  values  of  the  force. 
This  is  consistent  with  the  finding  of  Ladd  (1991).  In  the  range  of  5<r<7,  the  stress 
integration  method  gives  more  scattered  result  than  the  method  of  momentum  exchange. 
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For  smaller  radius,  i.e.,  coarser  lattice  resolution,  while  both  methods  give  poor  results 
(due  to  insufficient  resolution),  the  stress  integration  yields  much  larger  errors. 

Figure  4.6b  compares  Co  obtained  from  the  methods  of  momentum  exchange  and  the 
stress  integration  for  Re=10.  The  momentum  exchange  method  seems  to  gives  a 
converged  result  at  larger  r (>8).  Based  on  the  data  for  r>8,  an  average  values  of  CD 
-3.356  is  obtained.  In  contrast,  the  stress  integration  method  has  a larger  scattering  than 
the  large  r result  from  the  momentum  exchange  method  even  for  r>8.  Averaging  over 
the  result  for  r>8,  value,  the  stress  integration  gives  Co  —3.319.  The  difference  between 
converged  results  of  two  methods  is  about  1%.  For  r less  than  or  around  5,  the  scattering 
in  Co  from  the  stress  integration  method  is  much  larger  than  that  in  the  momentum 
exchange  method.  The  conclusions  from  the  comparisons  in  Figure  4.7  are  as  follows:  i) 
both  methods  for  force  evaluation  can  give  accurate  results;  ii)  the  momentum  exchange 
method  gives  more  consistent  drag;  iii)  in  the  range  of  10  < Re  < 100,  a resolution  of  10 
lattices  across  the  diameter  of  the  cylinder  are  needed  in  order  to  obtain  consistent  and 
reliable  drag  values.  In  other  words,  the  lattice  (grid)  Reynolds  number  Re*  should  be 
less  than  1 0 in  the  calculations. 

It  is  worth  noting  that  the  wall  shear  stress  in  the  channel  flow  obtained  using  the 
method  of  momentum  exchange  has  a relative  error  that  is  proportional  to  the  resolution 
across  the  channel.  For  a resolution  of  10  —20  lattices  across  the  diameter  considered 
here,  the  relative  error  in  the  drag  appears,  however,  smaller  than  in  the  channel  flow 
case.  At  Re=100,  with  r>10,  the  average  value  of  the  drag  using  the  method  of 
momentum  method  has  a 1.7%  relative  error  comparing  with  Fomberg’s  data  (1991).  If 
the  boundary  layer  thickness  is  estimated  roughly  to  be  3x2r/Re'/2  - 6,  there  are  only 
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about  6 lattices  across  the  boundary  layer  over  which  the  velocity  profile  changes 
substantially.  Based  on  the  insight  from  the  channel  flow  result,  it  is  possible  that  the 
deviatoric  shear  stresses  on  the  surface  of  the  cylinder  that  are  effectively  incorporated  in 
the  method  of  momentum  exchange  suffer  comparable  level  of  error  as  in  the  channel 
flow.  The  effective  error  cancellation  over  the  entire  surface  of  the  body  may  have 
contributed  to  the  good  convergence  behavior  in  the  drag  shown  in  Figure  4.6a  and 
Figure  4.6b. 

4.3.3  Flow  over  an  Asymmetrically  Placed  Circular  Cylinder  in  Channel  with 
Vortex  Shedding 

Schafer  and  Turek  (1996)  reported  some  benchmark  results  for  a laminar  flow  over  a 
circular  cylinder  asymmetrically  placed  inside  a channel.  The  ratio  of  the  center-to- 
upper-wall  is  hjr= 2.1  and  the  ratio  of  the  center-to-lower-wall  is  hJr= 2.0.  In  the  LBE 
computation,  r=12.8  is  used.  The  computational  domain  is  564x105.  To  preserve  the 
geometric  similarity,  A of  the  upper  boundary  is  set  to  0.76  and  A of  the  lower  wall  is  0.2. 
The  relaxation  time  is  z=  055.  The  centerline  velocity  is  0.095.  The  channel  inlet  has  a 
parabolic  velocity  and  is  specified  at  2.5  radii  upstream  of  the  cylinder  center;  this  results 
in  A=0.2  at  for  the  inlet  boundary.  A zeroth-order  extrapolation  for  fa  is  used  at  the  exit. 

The  Reynolds  number  based  on  the  average  inlet  velocity  U and  cylinder  diameter  is 
tfe=100. 

At  this  Reynolds  number,  the  flow  becomes  unsteady  and  periodic  vortex  shedding  is 

I Fy  I 

observed.  Figure  4.7a  and  Figure  4.7b  compare  the  lift  coefficient,  Cl=  , and  the 

pU  r 

drag  coefficient  Cd  (see  Eq.  (4.13))  with  the  benchmark  results  (Schafer  and  Turek, 
1996).  We  first  note  that  the  numerical  value  of  Strouhal  number  St  = 2r  /(UT)  is  0.300 
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in  which  T is  the  period  of  the  lift  curve.  This  agrees  very  well  with  the  range  of  value 
(0.2995-0.305)  given  in  Schafer  and  Turek  (1996).  We  note  that  there  is  very  small 
difference  in  Cl  between  the  momentum  exchange  method  and  the  surface  stress 
integration  method.  For  the  drag  coefficient  Co(t),  it  is  interesting  to  note  that  although 
there  is  about  0.25%  difference  between  the  results  given  by  momentum  exchange 
method  and  the  surface  stress  integration  method,  both  methods  of  force  evaluation  give 
two  peaks  in  the  Co(t)  curves.  Physically,  these  two  peaks  in  Co(t)  curve  correspond  to 
the  a weaker  vortex  and  a stronger  vortex  behind  the  cylinder.  The  difference  in  the 
strength  of  the  vortices  results  from  the  difference:  hjr= 2.1  and  hJr= 2.0.  There  is  no 
report  on  the  occurrence  of  these  two  peak  in  Schafer  and  Turek’s  work.  Instead,  a range 
of  the  Co,max  (from  3.22  to  3.24)  was  given.  The  values  of  these  two  peaks  are  well 
within  this  range.  The  presence  of  two  peaks  may  be  the  sources  of  the  apparently  wider 
range  of  Co.max  than  Cl, max.  A further  refined  computation  of  the  present  problem 
using  a multi-block  procedure  (D.  Yu  et  al.  2001)  with  r=40  in  the  fine  grid  region  yield 
nearly  the  same  results  for  Co(t)  and  Cl(F). 

4.3.4  Pressure  Driven  Flow  in  a Circular  Pipe 

The  steady  state  flow  field  was  obtained  (Mei  et  al,  1999)  using  D3Q19  model  with 

z=0.52.  Eq.  (4.6)  is  used  to  evaluate  the  force  on  the  boundary  points  along  the 
circumference  of  the  pipe  over  a distance  of  one  lattice  in  the  axial  direction.  The 
resulting  axial  force,  Fx,  is,  equivalently,  the  force  given  by  tw{270a5k)  in  which  tw  is 

the  wall  shear  stress.  For  a fully  developed  flow  inside  a circular  pipe,  the  exact  fluid 
shear  stress  at  the  pipe  wall  is  given  by  the  following  relation: 

C"(*2  r)=^-w2. 

dx 


(4.14) 
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We  examine  the  normalized  axial  force, 

r1  = Fxl{^-nr2)  (4.15) 

ax 

Figure  4.8  shows  the  normalized  coefficient  r|  over  a range  of  r:  3.5~  23.5.  Except  for 
r< 5,  rj  is  quite  close  to  1.  It  was  noticed  in  Mei  et  al.  (2000)  that  the  accuracy  of  LBE 
solution  for  the  pipe  flow  is  not  as  good  as  that  for  the  2-D  channel  flow  due  to  the 
distribution  of  values  of  zl’s  around  the  pipe.  The  accuracy  of  the  drag  is  dictated  by  the 
accuracy  of  the  flow  field  if  the  force  evaluation  method  is  exact.  For  the  pipe  flow,  the 
error  in  Fx  results  from  the  inaccuracy  in  the  flow  field  and  the  errors  in  the  force 
evaluation  scheme  based  on  momentum  exchange  (as  seen  in  the  previous  section  for  the 
2-D  channel  flow  case).  For  r> 5,  the  largest  error  in  Fx  is  about  3.5%  and  it  occurs  at 
r=  15.5.  Again,  there  is  no  systematic  error  in  Fx.  Given  the  complexity  of  the  boundary 
in  this  3-D  case,  the  results  shown  in  Figure  4.8  are  satisfactory  in  that  it  adds  further 
credence  to  the  momentum  exchange  method  for  force  evaluation. 

4.3.5  Steady  Uniform  Flow  over  a Sphere 

To  limit  the  computational  effort,  a finite  domain:  -H/2<y<H/2  and  - HI2<z<HI2  with 
FHr=  10  is  used  to  compute  the  flow  over  a sphere  of  radius  r (Figure  4.9).  Two  cases  are 
considered.  One  is  to  simulate  an  unbounded  flow  over  the  sphere.  The  other  is  to 
simulate  flow  over  a planer  array  of  sphere  (all  located  at  x=0)  with  the  center  of  the 
spheres  forming  square  lattices.  In  the  former  case,  the  boundary  conditions  at  jy=  1 
(y=H/2  corresponds  to  jy  =2)  for /a’s  are  given  by  the  following  linear  extrapolation, 

fjjz,  1 Jx)  = 2fcJJz,  2,  /x)  -fjjz,  3 ,jx)  (4. 1 6) 


The  velocity  atjy=2  is  set  as 
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u(j'z,  2,  jx)  = u(jz,  3,  j x).  (4. 1 7) 

Similar  treatment  is  applied  aiy=H/2  and  z=±HI 2.  In  the  latter  case,  symmetry  conditions 
are  posed  on/a’s  at  to  jy  =1  using  the  values  of  fa's  at  jy  =3  (see  Mei  et  al.  1999)  for  the 
2-D  case).  At  the  inlet,  a uniform  velocity  profile  is  imposed  at yx=1.5  (half  way  between 
the  first  and  second  lattices).  The  upstream  boundary  is  located  at  7.5  radii  to  the  left  of 
the  sphere  center  in  all  simulations. 

For  flow  over  a sphere,  the  drag  coefficient  is  often  expressed  as 


CD  = 


-Fr 


1 rr2 

— pU  7W 

2 


- = —</>■ 
2 Re 


or 


0 = 


~FX 

6 jr/iUr 


(4.18) 


where  ^accounts  for  the  non-Stokesian  effect  of  the  drag.  For  two  types  of  the  boundary 
conditions  at  (y=±H/ 2,  z=±H/2),  ^ym  denotes  the  non-Stokesian  correction  for  the  case 
where  the  symmetry  conditions  are  imposed  at  (y=±H/2,  z=±H/2)  and  $,nb  denotes  the 
results  for  the  case  where  the  extrapolation  for  /a’s  is  used  at  (y=±H/ 2,  z=±H/2 ) in  order 
to  simulate  the  unbounded  flow. 

Figure  4.10a  shows  the  non-Stokesian  coefficient  $inb  for  r=3.0,  3.2,  3.4,  3.6,  3.8,  4.0, 
5.1,  5.2,  5.4,  5.6  and  5.8  for  Hlr  =10  at  Re  = 10.  The  relaxation  time  is  r=0.7.  With  this 
range,  the  number  of  the  boundary  nodes  on  the  surface  of  the  sphere  increases  roughly 
by  a factor  of  (5.8/3)  ~3.74;  the  actual  counts  of  the  boundary  nodes  x*  gives  a ratio 
2370/546=4.35.  The  largest  difference  is  1.9%  between  r=3.0  and  r=3.2  which  have  the 
least  resolution  in  the  series  investigated.  For  a uniform  flow  over  an  unbounded  sphere, 
an  independent  computation  using  finite  difference  method  based  on  the  vorticity-stream 
function  formulation  with  high  resolution  gives  a drag  coefficient  ^-1.7986  at  Re=10. 
The  largest  difference  between  this  result  and  the  LBE  results  is  1.36%  at  r=  3.2.  If  the 
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LBE  data  for  the  drag  is  averaged  over  the  range  of  r,  one  obtains  0-1.8086  which 
differs  from  1.7986  by  0.54%.  Hence,  the  LBE  solutions  based  on  r=3.0-5.8  give  quite 
consistent  drag  force. 

Figure  4.10b  shows  the  non-Stokesian  correction  factor  ^ym  for  a uniform  flow  over  a 
planar  array  of  spheres  for  r=3.0-5.8  for  Hlr  =10  and  Re  = 10.  It  is  important  to  note  that 
with  the  improvement  of  the  surface  resolution  by  a factor  of  4.35,  there  is  little 
systematic  variation  in  0,ym(r).  The  largest  deviation  from  the  average  value,  ^ym.avc 
-1.963,  is  1.1%  at  r=5.0.  It  is  clear  that  the  LBE  solution  gives  reliable  fluid  dynamic 
force  on  a sphere  at  r~ 3.5  for  a moderate  value  of  Re.  The  set  of  data  for  ^ym  is 
inherently  more  consistent  than  that  for  (j\ mb  since  the  boundary  symmetry  condition  can 
be  exactly  specified  at  (y=±H/2,  z=±H/ 2)  while  the  extrapolation  conditions  given  by  Eqs. 
4.16-17  do  not  guarantees  the  free  stream  condition  at  (y=±H/ 2,  z=±H/2).  Yet,  both  $,nb 
and  0iym  exhibit  remarkable  self-consistency  from  coarse  to  not-so-coarse  resolutions. 

4.4  Conclusions 

Two  methods  for  evaluating  the  fluid  force  in  conjunction  with  the  method  of  lattice 
Boltzmann  equation  for  solving  fluid  flows  involving  curved  geometry  have  been 
examined.  The  momentum  exchange  method  is  very  simple  to  implement.  It  is  shown  in 
the  channel  flow  simulation  that  momentum  exchange  method  is  not  an  exact  method. 
The  error  in  the  wall  shear  stress  is  inversely  proportional  to  the  resolution.  In  2-D  and  3- 
D flows  over  a bluff  body,  it  can  give  accurate  drag  value  when  there  are  at  least  1 0 
lattices  across  the  body  at  Re-100.  The  method  of  integrating  the  stresses  on  the  surface 
of  the  body  gives  similar  result  when  there  is  sufficient  resolution  but  a much  larger 
uncertainty  exists  when  the  resolution  is  limited  in  comparison  with  the  method  of 
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momentum  exchange.  In  addition,  this  method  requires  considerably  more  efforts  in 
implementing  the  extrapolation  and  integration  on  the  body  surface  in  comparison  with 
the  method  of  momentum  exchange.  The  method  of  momentum  exchange  is  thus 
recommended  for  force  evaluation  on  curved  bodies. 


Table  4.1  Comparison  of  fluid  stresses  at  y=H  in  a pressure  driven 

2-D  channel  flow  with  dpldx  = - 1 O'6,  Ny  = 35  and  z-  0.6. 
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Figure  4.1  Layout  of  the  regularly  spaced  lattices  and  curved  wall  boundary. 


Figure  4.2  Variation  of  pressure  coefficient  on  the  surface  of  a circular  cylinder  at 
Re=40,  HZr=  10. The  result  is  obtained  using  7=0.6.  radius=6.6  lattice  unit. 
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Figure  4.3  Lattice  distribution  in  channel  flow  simulations  with  arbitrary  A. 


Figure  4.4  Ratio  of  the  wall  force  evaluated  using  Eq.  (4.9)  to  the  exact  value,  % = 

M~^~\y=n/  [-  — — H ] in  2-D  channel  flow  as  a function  of  rfor  zl=0.2,  1/3,  0.5, 
dy  2 dx 

and  0.7. 
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Figure  4.5  Dimensionless  wall  slip  velocity  in  the  2-D  channel  flow  as  a 
functionof  rfor  zl=0.2,  1/3,  0.5,  and  0.7. 
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Figure  4.6b 

Figure  4.6  Drag  coefficient  for  a uniform  flow  over  a column  of  cylinder  over 
arange  of  radius  r a)  Re=100  and  b)  Re=10. 
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t (lattice  unit) 
Figure  4.7a 


t (lattice  unit ) 
Figure  4.7b 


Figure  4.7  Comparison  of  the  lift  coefficient  (Figure  4.7a  ) and  drag  coefficient 
(Figure  4.7b)  with  the  benchmark  results  given  in  Schafer  and  Turek  (1996). 
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Figure  4.8  Ratio  of  the  tangential  force  Fx  on  the  pipe  to  the  exact  value  (^-7rr2) 

dx 


over  a range  of  pipe  radius  r. 


Figure  4.9  Schematic  for  uniform  flow  over  a sphere. 
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Figure  4.10a  Simulation  for  a single  sphere  in  an  unbounded  field 
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Figure  4.10b  Simulation  for  flow  over  a planer  array  of  spheres. 

Figure  4.10  Variation  of  the  non-Stokesian  correction  factor  with  sphere  radius  at 
Re=10. 


CHAPTER  5 

MULTI-TIME-RELAXATION  MODELS  IN  THE  LBE  METHOD 

5.1  Introduction 

In  Chapters  2 and  3,  the  single-relaxation-time  (SRT)  model  was  focused  on.  In 
attempting  to  obtain  solutions  for  high  Reynolds  number  flows  using  the  LBE  method,  It 
was  found  that  the  solution  field  for  (p,  u,  v)  often  exhibit  spatial  oscillations  in  regions  of 
large  gradients  such  as  stagnation  point  and  sharp  convex  corners.  Especially  near  a 
sharp  convex  corner,  because  the  pressure  and  the  vorticity  are  singular  locally,  a large 
gradient  in  the  density  or  pressure  field  exists.  Since  there  usually  is  insufficient 
resolution  near  the  corner,  the  large  gradient  is  often  accompanied  by  spatial  oscillations. 
Depending  on  the  geometry,  such  spatial  oscillation  can  propagate  into  the  flow  to 
contaminate  the  macroscopic  variables  in  a large  region  of  interest.  The  spatial 
oscillation  often  adversely  affect  the  computational  stability  and  convergence  rate. 

Recently,  Lallemand  and  Luo  (2000)  performed  detailed  analyses  on  the  dispersion, 
dissipation,  and  stability  characteristics  of  a generalized  lattice  Boltzmann  equation 
model  proposed  by  d’Humieres  (1992).  It  was  found  that  by  the  use  of  multiple 
relaxation  times  in  the  generalized  lattice  Boltzmann  equations,  better  computational 
stability  can  be  achieved  over  the  standard  LBGK  scheme  due  to  the  separation  of  the 
relaxations  of  the  various  kinetic  modes  in  the  generalized  lattice  Boltzmann  equation 
model— hereinafter  referred  as  the  multi-relaxation-time  (MRT)  model.  It  was  also  found 
in  Lallemand  and  Luo  (2000)  through  the  linearized  analysis  on  the  MRT  model  for 
various  simple  flows  that  the  MRT  model  gives  the  same  results,  to  the  second  order 
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accuracy,  as  the  SRT  LBGK  model  does.  It  seems  that  these  two  models  are  equivalent 
in  the  long  wavelength  limit  for  macroscopic  variable  of  interest  and  the  difference  is  a 
high  order  effect  based  on  their  linear  analysis.  Such  high  order  differences,  however, 
can  be  hardly  detected  in  simple  linear  flows. 

Many  fluid  flow  problems  possess  mathematical  singularities.  Since  a singularity 
often  affects  numerical  solutions  in  high  wavenumbers,  it  is  expected  that  the  results  of 
the  MRT  model  be  noticeably  different  from  that  of  the  SRT  model,  at  least  locally  near 
the  geometric  singularity.  For  convection-dominated  problems,  such  local  difference  in 
the  solution  behavior  can  also  lead  to  difference  in  the  solutions  over  a larger  scale.  It  is 
insightful  to  investigate  how  the  solution  based  on  the  MRT  model  behaves  in  such  flows 
in  comparison  with  the  standard  LBGK  model. 

In  this  chapter,  detailed  comparison  of  the  performance  of  these  two  LBE  models  for 
various  flows  with  geometric  and  flow  singularities  are  provided.  A brief  background  on 
the  MRT  model  will  be  described  first.  Computational  results  for  the  pressure,  viscous 
stresses,  vorticity,  and  velocity  in  regions  of  large  gradient  will  be  compared  between  the 
MRT  and  SRT  models  under  otherwise  identical  computational  and  physical  parameters 
for:  1)  Stokes  first  problem,  2)  steady  uniform  flow  over  a cascade  of  zero-thickness, 
finite  length  flat  plates,  3)  steady  uniform  flow  over  a cascade  of  finite-thickness,  semi- 
infinite  length  plates,  and  4)  steady  lid-driven  cavity  flow.  The  flow  in  the  Stokes  first 
problem  is  singular  at  t=0.  The  flow  in  a lid-driven  cavity  has  two  corners  at  the 
intersection  of  the  moving  and  stationary  walls  in  which  the  viscous  stresses  have  non- 
integrable  singularities.  The  flows  over  a plate  and  a step  have  singularities  in  the 
pressure  and  stresses,  but  they  are  weaker  than  in  the  lid-driven  cavity  flow.  These  flows 
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with  varying  degree  of  singularities  allow  for  a detailed  assessment  of  the  two  LBE 
models.  As  will  be  demonstrated,  the  MRT  model  performs  better  in  flows  involving 
large  gradients  than  the  SRT  LBGK  model.  In  all  the  simulations,  the  initial  density  is  set 
to  be  po=  1.  The  results  for  the  density  (and  thus  pressure)  are  presented  only  in  terms  of 
the  deviation  from  po  or  some  upstream  reference  value.  The  results  of  the  SRT  model 
are  obtained  by  running  the  same  MRT  code  with  si  =sj,  =ss  =Sj  = s&  = sg  = 1/r. 
Obviously,  the  relative  performance  between  different  LBE  models  depends  on  many 
factors  including  the  solution  characteristics,  types  of  variables  under  investigation,  and 
grid  resolution.  No  exhaustive  comparison  will  be  made.  Instead,  we  have  chosen 
reasonable  grid  size  in  all  cases  to  contrast  the  behavior  of  the  two  models. 

5.2  Multi-Relaxation-Time  (MRT)  Model 
In  d’Humieres  ( 1 992),  a new  set  of  variables  R = (p,  e,  a,  jx  ,qx,jy,qy,pxx,  pxy ) ‘ are 

introduced  and  R is  related  to  the  set  of  F = (/„ , /, , f2 , f3 , /4 » fs » ft » /? » fs ) ‘ as  follows, 


where  M_  is  the  9x9  matrix  transforming  F to  R.  In  the  vector  R , p is  the  fluid  material 

density,  e is  the  energy,  a is  related  to  the  square  of  the  energy,  jx  and  jy  are  the 
momentum  density  (or  mass  flux),  qx  and  qy  correspond  to  the  energy  flux,  and  p ^ and  p^ 
correspond  to  the  diagonal  and  off-diagonal  component  of  the  viscous  stress  tensor.  One 


81 


of  the  inherent  disadvantage  of  the  SRT  LBGK  model  is  that  all  variables  are  relaxed  in 
the  same  manner  as  given  by  Eq.  (1.10a).  In  lieu  of  Eq.  (1.10a),  the  collision  procedure 


of  the  MRT  model  is  carried  out  as  follows, 

e = e-s2(e-eeq),  (5.2a) 

£ = £ -s3(£  -£eq) , (5.2b) 

9x=<ix-ss  (qx-qexq),  (5.2  c) 

Vy  =^y-Sy(qy-qeyv),  (5.2  d) 

Pxx  = pxx  - h (Pxx  - p2  ) > (5 .2  e) 

Pxy  = Pxy  ~ (Ply  - Ply  ) (5.2  f) 

where  ~ denotes  the  post-collision  state  and  the  equilibrium  values  were  chosen  to  be 

eeq  =-2p  + 3(u2  + v2),  (5.3a) 

£eq  =p-3(u2  +v2),  (5.3b) 

qexq=-u,  (5.3  c) 


= -v, 


(5.3  d) 


Pe*=u2-v2,  (5.3  e) 

P7y  = UV  (5.3  f) 

Before  the  streaming  step,  Eq.  (1.10b)  is  solved.  One  needs  to  transform  the  post- 


collision values,  R , back  to  F by  using 


F = M~]  R 


(5.4) 


In  practice,  Eq.  (5.4)  can  be  combined  with  Eq.  (5.2)  to  obtain  a single  expression 
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F =F-M~'S(R-R ) 
where  5 is  the  diagonal  matrix: 


(5.5) 


S = diag(  0,  s2 , s3 ,0,  s5 ,0,  s1,sa,s9). 


(5.6) 


The  streaming  step  in  the  MRT  model  is  carried  out  exactly  in  the  same  manner  for  each 
component  as  in  the  standard  LBGK  model  based  on  Eq.  (1.10b). 

In  Lallemand  and  Luo  (2000),  it  was  shown  that  for  the  MRT  model  to  give  the  same 
shear  viscosity  as  given  by  Eq.  (1.1 1)  for  the  SRT  model,  one  needs  to  set 

£8  = sg  = 1/r.  (5.7) 

It  is  more  flexible  to  chose  the  rest  of  the  relaxation  parameters:  S2,  S3,  s$,  and  s 7.  In 
general,  these  four  parameters  can  be  chosen  to  be  slightly  larger  than  1 . In  our  study,  we 
set  S2  =S3  =ss  =57  =1.2  for  simplicity.  Very  little  difference  is  observed  in  the  flow  field 
if  a different  value,  say  1.1,  is  used  for  (s2,  S3,  s$,  Sj).  It  is  worth  commenting  here  that  by 
setting  S2  =53  =55  =57  = Jg  = Sg  = 1/r,  the  SRT  model  is  recovered. 

5.3  Results  and  Discussions 


5.3.1  Stokes  First  Problem 

For  an  infinitely  long  wall  to  move  impulsively  with  a velocity  U at  t=0+,  the  exact 
solution  for  the  wall  shear  stress  is  given  by 


(5.8) 


where  v is  the  dynamic  viscosity  of  the  fluid. 
LBE  solutions  for  the  wall  shear  stress, 


I _ exact 

T 

I xy,w 


Figure  5.1  shows  the  relative  error  of  the 


exact  I 
xy,w  I 


(5.9) 
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The  results  are  obtained  using  the  relaxation  time  7=1Ax=0.55.  Near  t=0,  there  is 
insufficient  enough  spatial  resolution  for  the  Stokes  layer  of  the  thickness  fvt . Hence  as 
illustrated  in  Figure  5.1  substantial  oscillations  are  present  near  t=0.  Clearly,  the  error  in 
the  MRT  model  is  smaller  than  the  SRT  error  for  t<l  00  when  the  near-wall  velocity 
gradient  is  large.  Eventually,  the  effect  of  singularity  diminishes  and  the  two  solutions 
approach  each  other. 

5.3.2  Flow  over  a Cascade  of  Zero-Thickness,  Finite  Length  Flat  Plate 

A schematic  of  the  flow  is  shown  in  the  insert  of  Figure  5.2.  the  symmetry  condition 
is  imposed  at  y=±  H.  The  plate  is  placed  half-way  between  two  grid  lines  so  that  the 
conventional  bounce-back  condition  can  be  used  to  update  the  wall  condition  for /a’s.  A 
zeroth-order  extrapolation  is  used  at  the  downstream  exit  plane  for  fa’s.  A constant 
uniform  velocity  condition  is  imposed  at  the  inlet,  x/L=  -2.  The  plate  length  is  40  in 
lattice  unit  (by  taking  &=1).  The  relaxation  time  controlling  the  shear  viscosity  is  set  to 
be  r=0.512  and  the  Reynolds  number  based  on  the  length  is  Re=UL/v  =1000.  The  free 
stream  velocity  is  thus  (7=0.1  and  H/L=2  so  that  there  80  lattices  from  the  plate  to  the 
symmetry  line. 

Figure  5.2  compares  the  density  deviation,  p- 1,  as  a function  of  y at  x/L  =0.0125, 
which  is  half  grid  away  from  the  leading  edge,  based  on  both  the  MRT  and  SRT  models 
under  otherwise  identical  conditions.  Due  to  the  singularity  in  the  flow  at  the  leading 
edge,  it  is  inevitable  to  have  large  gradients  in  pressure,  stresses,  and  vorticity  near  the 
leading  edge  at  high  Re.  When  there  is  insufficient  numerical  resolution,  an  unphysical 
spatial  oscillation  develops  near  the  leading  edge.  However,  the  MRT  model  is  seen  to 
be  much  more  effective  in  suppressing  the  spatial  oscillation  for  p or  p near  the  leading 


84 


edge.  Figure  5.3  compares  p- 1 as  a function  of y atx/Z=0.5125  under  the  same  condition. 
Surprisingly,  the  solution  based  on  the  SRT  model  still  possesses  a substantial  level  of 
spatial  oscillations  even  in  the  middle  of  the  plate  for  the  whole  cross-section  while  the 
solution  from  the  MRT  model  has  become  sufficiently  smooth.  Figure  5.4  shows  the 
viscous  normal  stress,  r**,  normalized  by  pU/L,  as  a function  of  y at  x/L  = 0.5125. 
Similar  level  of  oscillations  is  observed  in  the  SRT  based  solution.  In  this  work,  the 
viscous  stresses  are  obtained  using  the  non-equilibrium  part  of  the  distribution  function 
as, 

V(1~)Z  l/.fcM)  4«.  •«.*,)  (5-10) 

zr  a=l  ^ 

Hence,  no  finite  difference  is  employed  for  the  evaluation  of  the  viscous  stresses.  Figure 
5.5  compares  the  dimensionless  viscous  shear  stress,  z^,,  as  a function  of  y at  x/L  = 
0.5125.  Again,  the  oscillations  in  the  SRT  based  solution  are  noticeable  outside  the 
viscous  boundary  layer. 

To  develop  a further,  quantitative  understanding  of  the  performance  of  the  two 
models  for  flow  over  a flat  plate,  the  streamwise  variation  of  various  macroscopic 
quantities  near  the  plate  ylL  =0.0125,  which  is  half-grid  above  the  plate,  are  also 
examined.  Figure  5.6  shows  the  variation  of  the  pressure  coefficient 


r P~P« 

C P „ Tt2 


PoU1 1 2 


(5.11) 


at  y/L  =0.0125  as  a function  of  x for  solutions  based  on  these  two  models  where  px  is  the 
pressure  at  the  centerline  of  the  inlet.  It  is  noted  that  the  singularity  at  x=0  resulted  in 
oscillation  in  Cp  for  about  4-5  grid  points  after  the  leading  edge  in  the  MRT  model. 
However,  Cp  in  the  SRT-based  solution  continues  to  oscillate  across  the  entire  plate. 
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Figure  5.7  shows  variation  of  the  viscous  normal  stress,  t**,  normalized  by  f^U/L,  aiy/L 
=0.0125.  The  superiority  of  the  MRT  model  over  the  SRT  model  can  be  clearly  observed 
in  regions  before  and  after  the  leading  edge.  Figure  5.8  compares  the  dimensionless  wall 
vorticity,  du/dy,  normalized  by  UIL,  between  the  two  models.  Little  oscillation  is 
observed  for  the  MRT  based  solution  while  the  SRT  based  solution  continues  to  show 
oscillatory  behavior  up  to  x/L  =0.4,  which  cannot  be  considered  as  the  local  region  of  the 
leading  edge. 

5.3.3  Lid-Driven  Cavity  Flow 

The  insert  in  Figure  5.9  shows  the  coordinate  system  for  the  flow  inside  the  cavity. 
The  first  line  of  the  grid  in  the  fluid  region  is  at  a distance  Adx  from  the  wall.  In  this 
study,  zl=0.3  is  used.  The  boundary  treatment  is  based  on  that  given  in  Mei  et  al.  (2000) 
for  curved  geometries.  The  height  of  the  cavity  is  H/8x=64+2A.  With  x =0.52  and 
/£e=1000,  the  velocity  of  the  moving  wall  needs  to  be  £7=0.1032. 

The  velocity  field  is  discontinuous  at  the  two  comers  on  the  moving  wall.  Thus  the 
flow  singularity  is  stronger  than  in  the  previous  two  cases  where  the  velocity  is 
continuous  near  the  convex  corner.  Figure  5.9  compares  the  x-component  velocity  as  a 
function  of  y at  x/H= 0.00464  which  is  on  the  first  grid  away  from  the  left  vertical  wall. 
Oscillations  are  observed  in  both  SRT-based  and  MRT-based  velocity  profiles  due  to 
insufficient  resolution  for  the  singularity.  However,  the  oscillation  in  the  MRT  solution 
has  smaller  amplitude  and  is  limited  to  a region  of  5-6  grids.  The  oscillation  in  the  SRT 
solution  has  larger  amplitude  and  propagates  further  into  the  flow  field.  Figure  5.10 
shows  the  vertical  component  of  the  velocity  as  a function  of  y at  x/H=  0.00464.  Again, 
the  SRT  solution  has  a much  larger  amplitude  and  larger  region  of  the  oscillation.  Figure 
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5.11  compares  the  velocity  profiles  of  the  x-component  of  the  two  solutions  at  the 
centerline  (x/H=  0.5)  in  the  lower  half  of  the  cavity  with  a finite  difference  solution  based 
on  the  vorticity-stream-function  formulation.  It  is  worth  noting  that  the  MRT-based 
solution  is  noticeably  more  accurate  than  the  SRT-based  solution  even  in  regions  where 
one  considers  far  away  from  the  singularities. 

As  a final  comment,  by  taking  advantage  of  many  zero  elements  in  M_~'  S and 

recognizing  various  common  factors  in  the  expressions  for  the  vector  the 

algorithm  for  the  collision  step  in  the  MRT  model  can  be  coded  quite  efficiently.  For  the 
entire  computation  of  the  collision,  streaming,  and  the  evaluation  of  macroscopic 
variables,  the  code  for  the  MRT  model  takes  only  about  10%  more  CPU  time  per  time 
step  than  an  SRT  code  does.  Flowever,  this  extra  10%  work  is  greatly  compensated  by 
the  improved  convergence  of  the  MRT  model  in  suppressing  efficiently  the  transient 
oscillation  associated  with  the  high-frequency  pressure  (acoustic)  waves,  the  much 
improved  quality  of  the  results,  and  the  reduced  demand  for  higher  resolution. 

5.4  Implementation  of  MRT  Model  in  the  Multi-Block  Method 

In  Chapter  3,  a multi-block  method  has  been  developed  for  the  single-relaxation-time 
(SRT)  model  in  the  LBE  method.  The  scheme  for  information  transfer  was  derived  to 
ensure  the  continuity  of  mass,  momentum,  and  stress  across  the  interface  between  the 
coarse  and  fine  blocks. 

In  the  preceding  section,  it  has  been  shown  that  the  MRT  model  can  reduce  pressure 
oscillation  in  the  solution.  A recent  study  for  the  3-D  MRT  shows  that  it  can  improve  the 
stability  in  computation  (d’Humieres  et  al.  2002  ). 
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In  what  follows,  a scheme  for  data  exchange  is  presented  for  the  multi-block  MRT 
model  and  computational  assessments  are  given. 

5.4.1  Basics  of  the  Multi-Block  Strategy  for  MRT  Model 

In  Chapter  3,  a multi-block  method  was  developed  based  on  consistency  of  viscous 
coefficients  on  different  blocks,  the  mass,  momentum  conservation  across  the  interface, 
and  continuity  of  stress  at  the  interface.  Consistency  of  viscous  coefficients  requires  that 
the  three  relaxation  time  s2,  s8,  s9  need  to  be  rescaled  according  the  same  rule  (Eq.  (3.3)) 
in  Chapter  3.  They  are  the  parameters  determined  the  viscosity  of  the  fluid.  s2  is 
relaxation  time  for  the  energy  mode,  and  s3,  s5,  and  s7  are  relaxation  times  for  energy 
square,  and  energy  fluxes.  Because  s2  is  rescaled,  it  is  required  that  s3,  s5,  and  s7  are 
also  need  to  be  rescaled.  Thus  all  the  relaxation  times  are  need  to  be  rescaled.  The  mass 
and  momentum  conservation  require  f{neq)  to  be  kept  unchanged  across  the  interface  and 

stress  continuity  requires  that  f{amq)  to  be  rescaled  at  the  interface.  For  the  MRT  model 
to  be  incorporated  in  the  multi-block  method,  a similar  strategy  for  the  transferring  of 
information  at  the  interface  needs  to  be  developed.  However  there  is  no  explicit 
expression  for  f^eq) in  the  MRT  model.  The  values  of  (e,e,qx,qy,pxx,pxy){eq)  are  known 
and  relationship  of  these  variables  to  velocity  and  density  are  given  by  Eq.  (5.3).  To 
implement  the  multi-block  method  base  on  the  MRT  strategy,  the  relations  between  f„q) 

and  (p, e, e,jx,qx,  jy ,qy,pxx, P*y )('9)  must  be  obtained. 

Substituting  / = f^eq)  into  Eq.  (5.1),  one  obtains 
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Suppose  that  the  equilibrium  distribution  function  satisfy  the  following 


(5.12) 


faq)=  wa[f*\  eau+  -^-(ea-u)2  - uu]  (5.13) 

c2  2c4  2c2 

which  is  an  incompressible  version  for  /je,)  (He  and  Luo,  1 997c)  over  the  conventional 
version  given  by  Eq.  (1.4)  for  the  SRT  model. 

Substitution  of  Eq.  (5.12)  into  (5.13)  yields  the  same  relation  between 
(< e,£,qx,qy,pxx,pxy)(eq)  and  the  velocity  and  density  as  given  by  Eq.  (5.3).  Hence,  it  is 
found  that 


P(eq)  = P 
Jxq)  = U 

j r=v 


(5.14) 


From  Eq.  (5.3)  and  Eq.  (5.14),  it  can  be  seen  that  to  maintain  the  continuity  of  u,  v, 
and  p across  the  interface,  one  needs. 


(. P . Jx » < lx  Pxx » Pxy)(eq’C)  = ( A e,  e,  j x ,q  x , j y,q  y , pa , Pxy  (5.15) 


Substituting  Eq.  (5.2)  into  Eq.  (5.4),  one  obtains 


F = M~'R-  M~'  S(R  - Rieq) ) 


(5.16) 


where  ~ denotes  the  post-collision  state  of  the  F vector 
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Using  R = R(‘‘n  + R{neq) , it  is  seen  that 


F = M~'  R(eq)  + M~'  (/  - §)R(neq) 


(5.17) 


where 
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(5.18) 


In  the  vector  R{"eq) , density  and  momentum  flux  do  not  appear  in  the  relaxation  or 


collision  procedure  even  if  the  corresponding  relaxation  parameter  in  S is  not  zero.  Thus 
these  parameters  can  be  any  value  and  they  will  not  affect  the  computational  result. 
Stresses  are  determined  by 


rv=  0 - §-  )(P„  ~ P? ) = 0 - y )P?*  (5.19) 

«•„=  (l-y)(^  -p2”)  = (1  -f)pLm>  (5.20) 


To  maintain  the  shear  stress  continuity  across  the  interface,  in  the  2-D  case,  it  requires 
that 
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(5.21) 


(5.22) 
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Other  components  in  the  vector  (p,e,£,jx,qx,jy,qy,pxx,pxy)(neq)  remain  unchanged 


across  the  interface  except  for  pxx  and  p . Thus  the  non-equilibrium  terms  in  the 
coarse  grid  can  be  expressed  using  the  fine  grid  results  as: 

0 

-i(f) 


R 


(neq,c) 


m- 


ss 

So 


m- 


e-e(eq) 

o 

£-£(eq) 

e-eM 

0 

£-£{eq) 

0 

0 

Y’(f) 

9.-9?" 

9.-9“ 

0 

(/) 

8 (P  P^) 

(c)  yy xx  v xx  ) 

8 

9,-9?" 

(p„-  p'S') 

(/) 

9 (p  p (e9)) 

(c)  SPxy  Pxy  > 

_(P*y-P%q)) 

(5.23) 


where 


„(/)  „(/) 

TU)  = diag(\,\,l,\ ,1,1,1, 

58  s9 


Now  consideration  is  given  to  the  post-collision  state.  In  coarse  grid  block, 


F(c)  = M~l  R{eq'c)  +M~\l-  S(c) )R{neq-c) 


Using  the  condition  that  Rieq'c)  = R{eqJ  ’ and  Eq.  (5.23),  one  obtains 

F(c)  =M~'R(eq’f)  + M~\l-S(c))T(f)R(neq'f) 

The  non-equilibrium  part  of  the  distribution  function  in  fine  grid  is  given  by 


j^(neQ,f)  g(/K-l 


MFif)  -R(eqJ) 


Substituting  Eq.  (5.27)into  Eq.  (5.25)  yields 

F(c)  =M~]R{ev'f)  + M"'(7  - S{c))T(f)  (7  - 5(/)r‘ 


MF(f)  -R(eqJ) 


(5.24) 


(5.25) 


(5.26) 


(5.27) 


(5.28) 
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In  transferring  the  data  from  a coarse  grid  to  a fine  grid,  one  similarly  obtains 


Fif)  = M ' R(eq’c) 


+ M~\l-  S(f) )T(C)  (7  - S(c) MF(C)  - R(eq'c)  I (5.29) 


where 


(5.30) 


5.4.2  Computational  Assessment 

5.4.2. 1 Steady  flow  over  the  NACA0012  airfoil 

Flows  over  the  NACA  0012  airfoil  (Figure  3.20)  at  Re=5000  are  computed  with  the 
present  multi-block  LBE  scheme.  The  computational  domain  and  block  configuration  are 
the  same  as  the  ones  in  Chapter  3.  In  the  first  case,  s8=s9=l/0.502,  s2=1.81,  and 
s3=s5=s7=1.2  are  given  in  coarsest  blocks.  s2,  s8,  and  s9  are  rescaled  in  different 
blocks.  s3,  s5,  and  s7  are  kept  unchanged.  Figure  5.12  shows  the  result  of  drag 
coefficient.  It  can  see  that  the  computation  is  not  converged.  Lallemand  and  Luo  (2000) 
suggest  that  s3,  s5,  and  s 7 should  be  chosen  slightly  larger  than  1 to  reduced  the  over- 
relaxation effects  of  these  modes.  To  test  if  that  requirement  is  necessary  in  multi-block 
method,  a second  case  is  computed.  In  the  second  case,  s8=s9=l/0.502,  s2=1.81,  and 
s3=s5=s7=1.2  are  given  in  the  coarsest  grid.  All  the  relaxation  times  are  rescaled.  We 
have  s3=s5=s7=0.85  in  the  intermediate  grid  and  s3=s5=s7=0.32  in  the  finest  grid. 
Figure  5.16  shows  the  non-convergent  drag  coefficient.  In  the  final  case  s8=s9=l/0.502 , 
and  s2=s3=s5=s7=1.8  are  given  in  the  coarsest  grid.  All  the  relaxation  times  are 
rescaled.  s2,  s3,  sJ,  and  s7  are  equal  to  1.09  in  the  finest  block.  Figure  5.14  shows  the 
convergent  history  of  drag  coefficient.  In  this  case,  computation  converges  to  steady 
state.  The  final  value  of  drag  is  0.05379  comparing  with  the  value  of  0.05409  by  Drela 
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and  Giles  (1987).  It  can  be  concluded  from  above  cases  that  all  relaxation  times  are  need 
to  be  rescaled  and  s3,s5,  and  s7  are  need  to  be  larger  than  1 when  different  values  are 
assigned  to  different  relaxation  times. 

5.4.2. 2 Lid-driven  cavity  flow 

The  lid-drive  cavity  flow  is  used  to  assess  the  present  multi-block  strategy  with  the 
MRT  model,  especially  for  the  interface  conditions.  The  computational  domain  and 
parameter  are  identical  to  the  ones  used  in  Chapter  3.  The  relaxation  parameter  are  S2  =sj 
=Ss  =sj  =1.8.  The  position  of  the  center  of  the  primary  vortices  is  (0.6172,  0.7350)  for 
Re=100,  compared  well  with  the  value  (0.6172,  0.7344)  from  Ghia  et  al.  (1982).  The  u- 
and  v-components  of  the  velocity  along  the  vertical  line  and  horizontal  line  through  the 
geometry  center  are  shown  in  Figure  5.15a  and  Figure  5.15b,  respectively.  It  is  seen  that 
the  velocity  profile  by  both  MRT  and  SRT  models  agree  well  with  results  from  Ghia  et 
al.  (1982).  Figure  5.16  shows  the  pressure  contours  obtained  from  the  multi-block 
solution.  Figure  5.17  shows  a local,  enlarged  view  of  the  pressure  contour  around  an 
interface  corner  point  indicated  by  the  circle  in  Figure  5.16.  Clearly,  the  pressure  is 
smooth  across  the  interface.  Figure  5.18-Figure  5.20  show  the  contours  of  shear  stress 
and  velocity  components,  which  are  also  smooth  across  the  interface. 

To  demonstrate  more  clearly  the  conservation  and  continuity  characteristics  at  the 
interface,  the  velocity,  density,  and  shear  stress  along  the  interface  AB  in  Figure  5.16  are 
shown.  Again,  the  values  of  physical  variables  from  the  fine  block  on  interface  AB  are 
obtained  by  second  order  extrapolation  from  inside  the  fine  grid  domain.  Figure  5.21  and 
Figure  5.22  show  the  velocity  on  the  interface.  It  can  be  seen  that  the  velocity  profiles 
between  coarse  and  fine  grids  agree  very  well.  Figure  5.23  shows  pressure  on  the 
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interface  which  are  smoother  in  both  coarse  and  fine  blocks  than  those  suing  the  SRT 
model,  shown  in  Figure  3.17.  However  there  is  a larger  discrepancy  of  pressure  on 
interface  between  fine  and  coarse  grid  near  the  upper  wall.  It  is  believed  that  the  source 
of  this  discrepancy  is  from  the  treatment  of  the  upper  wall  boundary.  As  mentioned  in 
Chapter  3,  the  different  A at  the  intersection  of  two  blocks  at  upper  wall  may  cause 
inaccuracy  in  the  boundary  treatment.  It  is  also  interesting  to  note  that  the  oscillation  of 
pressure  in  coarse  grid  did  not  transfer  into  fine  grid. 

5.5  Conclusions 

Based  on  the  detailed  examination  of  the  flow  fields  in  selected  cases,  it  is  clear  that  the 
MRT  model  has  substantial  advantages  over  the  SRT  model  in  handling  the  geometric 
singularities.  The  MRT  model  in  general  provides  smoother  variations  of  the 
macroscopic  quantities  and  has  much  smaller  regions  of  the  oscillation  near  a singularity. 
Since  the  spatial  oscillation  is  often  accompanied  by  the  high  frequency  pressure  waves 
in  transient  simulations,  the  MRT  model  also  offers  a better  convergence  toward  steady 
state  as  well.  The  MRT  model  is  strongly  recommended. 

A multi-block  scheme  is  developed  for  the  MRT  model  in  the  LBE  method.  The 
interface  condition  for  the  MRT  model  is  derived  to  ensure  the  mass  conservation  and 
stress  continuity  between  neighboring  blocks.  This  new  treatment  retains  the  inherent 
advantages  of  both  MRT  model  and  the  multi-block  methods. 
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Figure  5.1  Comparison  of  the  relative  error  in  the  evolution  of  the  wall  shearstress 
for  Stokes  first  problem  between  the  SRT  and  MRT  models. 


Figure  5.2  Comparison  of  the  density  profiles  near  the  leading  edge  (x/L=0.0125) 
between  the  SRT  model  and  MRT  model  at  Re=1000. 
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Figure  5.3  Comparison  of  the  density  profiles  at  x/L=0.5125  between  the  SRT 
model  and  MRT  model  at  Re=1000. 


Figure  5.4  Comparison  of  the  viscous  normal  stress  profiles  at  x/L=0.5125 
between  the  SRT  model  and  MRT  model  at  Re=1000. 
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Figure  5.5  Comparison  of  the  viscous  shear  stress  profiles  at  x/L- 0.5125  between 
the  SRT  model  and  MRT  model  at  /?e=1000. 


Figure  5.6  Comparison  of  pressure  coefficient  as  a function  of  x at  _y/Z.=0.0125 
between  the  SRT  model  and  MRT  model. 
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Figure  5.7  Comparison  of  viscous  shear  stress  as  a function  of  x at  y/L= 0.0125 
between  the  SRT  model  and  MRT  model. 
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Figure  5.8  Comparison  of  wall  vorticity  as  a function  of  x between  the  SRT  model 
and  MRT  model. 
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Figure  5.9  Comparison  of  the  velocity  profiles  of  x-component  at  x=A  (7=2). 


Figure  5.10  Comparison  of  the  velocity  profiles  of  y-component  at  x=A  (i= 2). 


99 


Figure  5.1 1 Comparison  of  the  velocity  profiles  of  x-component  at  x/H=  0.5  in  the 
lower  region  of  the  cavity. 


Figure  5.12  Convergent  history  of  drag  coefficient,  here  s8=s9=l/0.502 , s2=1.81, 
s3=s5=s7=1.2  in  coarsest  block.  s3,  s5,  and  s7  are  not  rescaled. 
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Figure  5.13  Convergent  history  of  drag  coefficient,  here  s8=s9=l/0.502,  s2=1.81, 
s3=s5=s7=l  .2  in  coarsest  block.  All  parameters  are  not  rescaled. 


Figure  5.14  Convergent  history  of  drag  coefficient,  here  s8=s9=l/0.502, 
s2=s3=s5=s7=1.81  in  coarsest  block.  All  parameters  are  not  rescaled. 
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Figure  5.15a  Comparison  of  u-velocity  along  the  vertical  line  through  geometric 
center. 


Figure  5.15b  Comparison  of  v- velocity  along  the  horizontal  line  through 
geometric  center. 

Figure  5.15  Comparison  of  velocity  between  SRT  model,  MRT  model  and  those 
by  Ghia  et  al.  (1982). 
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Figure  5.16  Pressure  contours  in  the  cavity  from  the  multi-block  MRT  LBE 
solution.  (For  the  circled  region,  see  Figure  5.17) 


Figure  5.17  Enlarged  view  of  pressure  contour  in  the  circled  region  in  Figure 
5.16  near  the  intersection  of  three  blocks.  The  figure  demonstrates  that  the 
block  interface  and  corner  are  well  handled  in  MRT  model. 
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Figure  5.18  Shear  stress  contour  based  on  the  MRT  model.  Solid  and  dash 
lines  represent  positive  and  negtive  values,  respectrively. 


Figure  5.19  Contour  of  x-component  of  velocity  based  on  the  MRT  model. 
Solid  and  dash  lines  represent  positive  and  negtive  values,  respectrively 
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Figure  5.20  Contour  of  y-component  of  velocity  based  on  the  MRT  model. 
Solid  and  dash  lines  represent  positive  and  negtive  values,  respectrively 


Figure  5.21  The  x-component  of  the  velocity  ux/U  on  the  interface  AB  based  on 
the  MRT  model.  In  Figure  5.21-Figure  5.22,  LM0.0156. 
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Figure  5.22  The  x-component  of  the  velocity  uy  /U  on  the  interface  AB  based  on 
the  MRT  model. 


Figure  5.23  Pressure  on  the  interface  AB  based  on  the  MRT  model. 
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Figure  5.24  Shear  stress  Txy  /(juU/H) on  the  interface  AB  based  on  the  MRT  model. 


CHAPTER  6 

IMPROVEMENT  OF  BOUNDARY  TREAMENTS  IN  THE  LBE  METHOD 
Two  classes  of  boundaries  in  LBE  simulation  are  considered  here:  the  open  boundary 
and  solid  wall  boundary.  Most  of  the  previous  research  efforts  have  focused  on  the 
treatment  of  the  solid  wall  boundary.  Here  we  first  give  a brief  description  on  the 
development  of  solid  wall  boundary  conditions. 

The  bounce-back  condition  originated  from  lattice  gas  automata  has  been  extensively 
used  in  LBE  simulation  (Ziegler  1993,  Ginzbourg  and  Alder  1994).  The  bounce-back 
scheme  in  the  LBE  method  is  intuitively  derived  from  LGA.  In  this  scheme,  the  particle 
(the  basic  variable  in  LGA)  hits  the  wall  and  bounces  back  opposite  to  its  incoming 
direction.  In  LBE,  the  particle  number  is  replaced  with  a real  variable  fa.  This  scheme  is 
very  simple,  efficient  and  easy  to  implement.  On  the  other  hand,  the  bounce-back 
scheme  is  only  a first  order  treatment  in  general  (Ginzbourg  and  Alder  1994).  In 
simulating  suspension  flow,  Ladd  (1994)  placed  the  wall  halfway  between  the  nodes  and 
it  has  been  shown  that  this  treatment  has  the  second  order  accuracy  for  the  straight  wall. 

The  bounce-back  scheme  has  been  modified  to  reduce  the  wall  slip  velocity  by  Nobel 
et  al.  (1995),  Inamuro  et  al.  (1995),  and  Maier  et  al.  (1996).  Zou  and  He  (1997)  extend 
the  bounce-back  condition  to  non-equilibrium  part  of  distribution  function. 

To  improve  the  numerical  accuracy,  Skordos  (1992)  suggest  adding  a first  order  non- 
equilibrium distribution  in  the  equilibrium  distribution  at  the  wall  nodes.  Chen  et  al. 
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(1996)  placed  the  wall  at  the  node  and  used  extrapolation  of  fa  on  the  fluid  side  to  obtain 
the  fa  at  x*. 

All  the  boundary  treatments  mentioned  above  modeled  curved  wall  as  zig-zag  steps, 
which  results  in  geometric  discontinuities  and  affect  the  computational  accuracy.  This 
error  is  amplified  when  Re  increases. 

Filippova  and  Hand  (1998)  presented  a curved  boundary  condition  using  Taylor 
series  expansion  in  both  space  and  time  for  fa  near  the  wall.  This  boundary  condition 
provides  a second  order  accurate  treatment  for  a curved  solid  wall  boundary.  In  addition, 
two  other  boundary  treatments,  one  by  Mei  et  al.  (1999)  and  the  other  by  Bouzidi  et  al. 
(2001)  have  also  been  presented.  The  details  of  above  three  boundary  conditions  will  be 
provided  in  Section  6.2  along  with  a new  solid  wall  boundary  treatment. 

In  general,  open  boundaries  refer  to  lines  (or  planes)  of  symmetry,  periodic  boundary, 
infinity,  and  inlet  and  outlet.  For  symmetric  and  periodic  boundaries,  the  condition  can  be 
given  exactly.  There  is  not  much  attention  given  to  the  impact  of  the  inlet  condition  on 
the  quality  of  the  solution  in  LBE.  Usually  the  solid  wall  boundary  condition  is  used  to 
provide  the  distribution  functions  at  inlet.  This  treatment  can  cause  instability  in  the 
computation,  because  the  fluctuation  inside  the  computational  domain  will  be  reflected 
back  from  the  inlet.  In  Section  6.1,  a new  inlet  treatment  is  presented  to  address  this 
issue. 

When  a boundary  moves,  new  issues  arise.  For  example  a fluid  node  can  experience 
the  change  of  phase.  These  issues  will  be  discussed  in  Section  6.3  with  presenting  a 
moving  wall  boundary  condition. 


109 


6.1  The  Open  Boundary  Treatment 

6.1.1  Introduction 

Different  methods  have  been  implemented  to  obtain  the  distribution  function  fa  at 

the  inlet,  from  which  the  specified  velocity  or  pressure  can  be  recovered.  The  most 
commonly  used  approach  is  based  on  the  bounce-back  of  the  distribution  function  or  non- 
equilibrium part  of  distribution  function  (Cornubert  et  al.  1991,  Ziegler  1993,  Behrend 
1995).  Typically  one  places  the  inlet  boundary  half  way  between  the  inlet  boundary  grid 
and  the  grid  after  the  inlet  with  a known  velocity  profile  specified  at  the  inlet.  The 
standard  bounce-back  scheme  for  fa  at  the  inlet  gives: 

~ 3 

fa  inlet  ~ fa  + 2\V ap  i„temal  2 inlet  (6-1) 

C 

where  wa  is  the  weighting  factor,  and  e-  = -ea  . 

Grunau  (1993)  placed  the  inlet  at  the  lattice  node  as  opposed  to  the  lattice  link  half- 
way between  the  grids.  He  then  assigned  equilibrium  distribution  function  to  be  the 
desired/at  the  inlet.  For  the  2-D,  9-velocity  model,  this  method  gives: 

3 9 3 

fa  inlet  =faq)  = Pinlet  [1  + — , • «,„/*  + ~ T (*a  ' Y ~ ~Y  U*  inlet]  (6.2) 

c 2 c 2c 

This  method  was  considered  to  have  a larger  error  than  the  one  using  Eq.  (6.1)  (Zou  and 
He  1997). 

Skordos  (1992)  presented  a scheme  to  calculate  fa  from  a given  initial  velocity  and 
density  field  (1992).  In  this  scheme  the  boundary  is  placed  at  the  lattice  nodes  and  the 
first  order  term  in  the  Chapman-Enskog  expansion  /(1),  is  added  to  f(eq),  to  obtain  the 


desired  distribution  function, 


In  computing  /(l> , the  gradients  of  velocity  and  density  are  required.  For  computational 
purpose,  the  gradient  was  obtained  by  using  finite  difference  inside  the  computational 
domain  near  the  inlet.  If  the  inlet  velocity  derivative  is  known,  there  is  no  need  to  use 
finite  difference  to  obtain  /(l) . This  method  can  result  in  very  stable  computation  if  the 
analytical  value,  instead  of  the  finite  difference,  for  the  gradient  is  used.  When  finite 
difference  is  used  in  constructing  /(l) , since  the  value  inside  the  computational  domain  is 
used  computational  stability  may  be  compromised. 

In  Eq.  (6.1),  fa  at  the  inlet  is  determined  using  the  variables  of  f’s  inside  the 
computational  domain.  Hence, 

fa_mie,  = Function(f a interml , pMemal , umlel ) (6.4) 

In  Eq.  (6.3)  the  equilibrium  part  of  distribution  function  is  determined  by  the  values  of 
known  density  and  velocity  at  the  inlet, 

fsZ lie,  = Fu”ctior{p,nlel , Uinlel ) (6.5) 

If  the  gradients  for  the  density  and  velocity  are  evaluated  using  finite  difference,  the  non- 
equilibrium part  is  related  to  the  values  the  velocity  and  density  inside  the  computational 
domain  so  that: 

/i%  = FunCti°n(P internal  > “.ntema,  ) (6-6) 

If  the  gradients  are  given  analytically,  fF)mk,  does  not  vary  with  the  flow  variables 

inside  the  domain.  However,  this  can  result  in  over-specification  of  the  boundary 
conditions  (for  macroscopic  variables)  as  one  needs  both  the  value  and  derivative  of  the 


Ill 


variable  on  one  side  of  the  domain.  When  the  inlet  velocity  is  uniform,  Eq.  (6.3)  can  be 
further  simplified  to 


f _ f(eq) 

J a inlet  J a inlet 


(6.7) 


We  refer  to  Eq.  (6.7)  as  the  equilibrium  inlet  condition. 

We  now  can  see  that  equations  (6.1)  and  (6.3)  are  different  in  one  major  aspect.  In 
Eq.  (6.1),  f-  mlel  is  the  function  of  p and  fa  ’s  of  the  interior  nodes  whereas  in  Eq.  (6.4), 


f-  miei  's  onb'  function  of  gradient  of  the  p and  fa  ’s. 

In  Skordos’s  approach,  the  specification  of  the  gradients  of  density  and  velocity 
analytically  eliminates  the  interaction/coupling  between  the  boundary  nodes  and  the 
interior  nodes,  but  the  coupling  can  get  stronger  when  the  finite  difference  is  used  to 
evaluate  the  derivatives  since  the  values  at  the  interior  points  are  used.  Such  a strong 
coupling  often  results  in  strong  pressure  wave  to  originate  from  the  inlet. 

In  most  calculations  using  Eq.  (6.1)  to  specify  the  inlet  condition,  the  pressure  field 
exhibits  fluctuation  during  calculation.  This  can  be  understood  now  since  the  bounce- 
back  condition  given  by  Eq.  (6.1)  allows  for  strong  coupling/interaction  between  the 
interior  and  the  inlet  boundary.  This  kind  of  interaction/coupling  will  become  a serious 
problem  when  computing  high  Reynolds  number  flow.  Figure  6.1  shows  the  pressure 
contour  for  a flow  over  NACA0012  airfoil  at  Re=2000.  The  boundary  conditions  at  inlet, 
upper  and  lower  boundaries  are  given  using  bounce-back  condition.  It  can  easy  see  that 
there  is  a strong  interaction  between  inlet  and  interior  field.  This  interaction  is  intensified 
at  the  two  corners  near  the  inlet,  because  at  this  region  the  interaction  comes  from  both 
the  inlet  and  the  upper  or  lower  boundary. 
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6.1.2  The  Proposed  Inlet  Treatment 

Equation  (6.3)  gives  an  inlet  condition,  which  does  not  reflect  the  pressure  oscillation 
back  into  computational  domain,  but  the  finite  difference  used  to  obtain  /(l)  can  cause 
computational  instability  (Skordos,  1992).  Here  we  propose  an  inlet  boundary  condition 
that  retains  the  basic  characteristic  of  Eq.  (6.3);  that  is  the  equilibrium  part  of  the 
distribution  function  at  inlet  is  obtained  using  the  known  inlet  velocity  and  density,  while 
/(l)  is  replaced  with  non-equilibrium  part  of  distribution  function  f(noneq) : 

f _ f(eq)  , r(nomq) 

J a _ inlet  J a _ inlet  ' J a _ inlet  \0. 

Here  a 1-D  example  is  used  to  clarify  the  how  to  find  the  values  of  f^eq]nlet  and 
in  Eq.  (6.8),  which  correspond  the  fxeqB  and  f"°neq  in  Figure  6.3  . Considering 
post  streaming  stage  at  the  inlet,  boundary  conditions  need  to  be  given  to  specified  /,  B , 
which  include  two  components:  fxqB  and  fx°Beq . At  the  inlet  fluid  density  and  velocity 
are  known,  thus  the  equilibrium  distribution  function  fxeq)  at  inlet  can  be  determined. 
Knowing  the  fxq  and  fxeq , fxqB  can  be  obtained  easily  using  liner  interpolation: 


f(«7)  _ , r(eq)  _ 

J\,B  ~ J\J  J],l  ) , . 

1 + A 


The  inlet  can  be  placed  arbitrary.  It  takes  a little  bit  more  effort  to  find 


using  an  extrapolation  to  find  f"'jncq  , such  as 


(6.9) 

fxJeq . To  avoid 


/•  noneq  r noneq  _|_  / r noneq  r noneq  \ a 

5 ,/  J 5 ,B  5 ,B  J 5 ,C  * 

which  maybe  lead  the  instability  in  computation,  we  let 

/*  noneq  r noneq 

5,1  ~ J5,B 


(6.10) 


(6.11) 


113 


Then  momentum  balance  requires 


r noneq 


r noneq 
1,5 


(6.12) 


J\j  i 


Finally  an  interpolation  is  used  to  obtain  /, 


r noneq  , 
1 ,B 


r noneq 
\,B 


r noneq 
1 ,C 


r noneq  \ 

Ju  )-T+a 


noneq  \ 


(6.13) 


Using  this  method,  the  non-equilibrium  part  of  the  distribution  function  can  be  added 
to  the  equilibrium  distribution  function.  We  call  this  inlet  boundary  condition  as  local 
inlet  condition 

6.1.3  Computational  Assessment 

In  the  computational  cases  to  be  presented  and  discussed,  the  following  questions  are 
addressed: 

1.  Will  the  local  inlet  condition  give  the  same  accurate  results  as  those  of  bounce-back  if 
both  conditions  lead  to  converged  results? 

2.  How  does  the  inlet/free-stream  boundary  condition  impact:  the  convergence  of  the 
solution  and  the  quality  of  the  flow  field?  Flows  over  a column  of  circular  cylinders 
under  different  initial  condition  and  different  Reynolds  numbers  and  a flow  over  an 
airfoil  are  used  to  assess  these  issues. 

6. 1.3.1  Steady  uniform  flow  over  a column  of  cylinders  under  /?e=40 

For  a uniform  flow  over  a column  of  circular  cylinders  of  radius  r with  center-to- 
center  distance  denoted  by  H,  symmetry  conditions  for/a’s  are  imposed  at  y=  ±H/2.  At 
the  inlet,  the  uniform  velocity,  u=V,  is  specified  using  bounce-back  boundary  condition 
and  the  local  inlet  condition. 


At  the  exit,  a simple  extrapolation  is  used, 
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UNx,j)  = 2UNx-\J)-UNx-2,j)  for  a=4,  5,  and  6.  (6.14) 

On  the  surface  of  the  circular  cylinder,  FH’s  boundary  condition  is  used  to  update  the 
fa  s;  see  Figure  6.2  for  the  computational  domain.  In  the  computation  z=0.52, 
Uiniet= 0.0127.  The  drag  coefficient  is  calculated  using  the  momentum  exchange  method. 

Figure  6.4a  shows  the  convergence  histories  of  the  drag  coefficient  using  these  two 
inlet  conditions.  It  is  seen  that  the  local  inlet  condition  results  in  much  smaller  oscillation 
and  faster  convergence  towards  the  steady  state  than  that  using  bounce-back  inlet 
condition.  Furthermore,  Cd  based  on  the  local  inlet  condition  possesses  much  smaller 
oscillation  than  that  using  bounce-back  inlet  condition  when  the  solution  is  near 
convergence.  This  means  that  the  pressure  oscillation  in  the  flow  field  is  much  small 
when  local  condition  is  used.  The  average  values  of  two  converged  results  are  nearly  the 
same  with  Cd  ~ 1 .6485  (Figure  6.4b). 

6.1.3. 2 Steady  uniform  flow  over  a column  of  cylinders  under  /?e=100  with  initial 
velocity  equal  inlet  velocity 

In  this  case  all  other  parameters  are  kept  the  same  as  those  in  the  first  case,  except 
that  Reynolds  number  is  increased  to  100  by  increasing  the  inlet  velocity.  Figure  6.5 
shows  the  convergence  results  after  vortex  beginning  to  shed.  The  bounce-back  condition 
gives: 

G/_average  =1.41 1 

C/_max=0.335 

The  local  inlet  condition  gives:  Q _average=l  -407 

C/_max=0.334 
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We  obtain  St= 0.33  for  both  cases.  This  case  shows  that  both  bounce-back  and  local  inlet 
condition  can  produce  same  results  if  convergence  is  achieved. 

6.1.3.3  Steady  uniform  flow  over  a column  of  cylinders  under  /fe=100  with  initial 

velocity  =0 

All  computational  parameters  are  the  same  as  those  in  case  2,  except  that  the  initial 
velocity  is  set  to  0,  which  means  that  there  will  be  a strong  interaction  between  inlet  and 
the  interior  of  the  flow  field.  Figure  6.6  shows  drag  and  lift  coefficients  using  inlet 
condition  based  on  the  bounce-back  condition.  Irregular  pattern  appears  in  C/  curve  after 
vortex  shedding  starts.  It  does  not  die  out  with  time.  It  suggests  that  strong  interaction 
exists  between  inlet  boundary  and  interior  field.  Figure  6.7  shows  results  using  local  inlet 
condition.  The  lift  coefficient  C/  behaves  quite  well.  It  suggests  a much  weaker 
interaction  between  inlet  boundary  and  interior  field. 

6.1.3. 4 A steady  flow  over  NACA0012  airfoil 

The  NACA  0012  airfoil  is  a popular  wing  model,  which  has  been  used  extensively. 
Flow  fields  at  ./?e=500,  1000,  2000,  and  5000  have  been  computed  using  the  multi-block 
LBE  scheme.  Details  have  been  given  in  Chapter  3.  At  the  inlet,  upper  and  lower 
boundaries,  local  inlet  condition  are  used  for  f's  based  on  the  free-stream  velocity,  At 
the  downstream  boundary  a zeroth  order  extrapolation  for  f's  is  used. 

The  bounce-back  inlet  condition  yields  unphysical,  large  scale  spatial-temporal 
oscillation  in  the  flow  field  when  Re= 2000  and  Re= 5000  ( 7=0.5025  and  0.50125 
respectively).  Figure  6.8  shows  the  variation  of  drag  with  time  under  Re=2000  using 
bounce-back  condition.  It  can  easily  see  that  flow  soon  become  unsteady  after  the 
computation  starts.  Figure  6.9  shows  the  variation  of  drag  with  time  under  Re=2000 
using  local  inlet  condition.  Using  local  inlet  condition  computation  convergences  very 
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faster  to  the  steady  state.  Cd=0.8601  is  very  close  to  the  value  of  0.8373  by  Xfoil.  Figure 
6.10  shows  the  density  contour  of  flow  field  using  bounce-back  at  time  equaling  to  30000 
lattice  unit.  The  irregular  and  large  scale  density  oscillation  shows  that  flow  field  has 
become  unphysical.  Figure  6.1 1 shows  density  contour  using  local  inlet  condition  at  time 
equaling  to  30000  lattice  unit.  It  can  see  that  there  is  no  oscillation  of  density  in  the  flow 
field.  Thus  the  local  inlet  condition  is  demonstrated  to  be  superior  over  the  bounce-back 
type  of  inlet  condition.  It  greatly  improved  the  computational  stability  and  the  quality  of 
the  solution. 

6.1.3. 5 Channel  flow  over  an  asymmetrically  placed  cylinder  at  Re=20 

Schafer  and  Turek  (1996)  reported  a study  of  a laminar  flow  over  a circular  cylinder 
placed  asymmetrically  inside  a channel.  In  the  present  study,  the  computational  domain 
and  parameters  are  the  same  as  the  ones  in  Section  4.3.3,  except  that  Re  equals  to  20 
present  case.  At  inlet  bounce-back  and  local  inlet  conditions  are  used.  Both  inlet 
conditions  give  Cd=5.59  which  is  in  the  range  of  5.57-5.59  by  Schafer  and  Turek  (1996). 
Figure  6.12  shows  the  comparison  of  velocity  profiles  along  the  y-direction  for  both  inlet 
conditions  at  nodes,  which  are  0.2  lattices  to  the  right  of  inlet.  Figure  6.13  shows  shear 
stress  profiles  along  the  same  direction  at  same  position  as  Figure  6.12.  It  can  see  that 
local  inlet  condition  can  successfully  specifies  the  stress  at  the  inlet.  In  Figure  6.13,  shear 
stress  profiles  are  not  linear  distributions  along  y-direction  because  of  the  affect  of 
cylinder. 

6.1.4  Conclusions 

Computations  based  on  the  bounce-back  inlet  condition  and  equilibrium  inlet 
conditions  are  performed  for  various  cases.  Bounce-back  inlet  condition  results  in  slower 
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convergence  toward  steady  or  dynamically  steady  state  and  causes  stronger  interaction 
between  the  inlet  and  interior  field.  Such  interaction  may  affect  the  long-term  behavior 
of  the  solution,  especially  at  higher  Re.  Using  the  local  equilibrium  value  at  the  inlet  can 
reduce  the  impact  of  the  boundary-interior  interaction  on  the  unsteady  development  of  the 
flow  field,  improves  the  convergence  rate,  improves  the  computational  stability,  and 
improves  quality  of  the  overall  solution. 


Figure  6.1  Density  contour  at  Re=2000  using  bounce-back  condition  at  inlet,  upper  and 

lower  boundaries 
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Figure  6.2  Computational  domain  and  boundary  conditions. 
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Figure  6.3b  Configuration  of  equilibrium  distribution  at  inlet. 
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Figure  6.3a  Configuration  of  non-equilibrium  distribution  at  inlet. 

Figure  6.3  Configuration  of  distribution  functions  which  are  used  to  construct  inlet 

boundary  condition. 
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Figure  6.4a  Convergence  history  of  the  drag  coefficient. 


Figure  6.4b  Drag  coefficient  near  convergent  state. 

Figure  6.4  Convergence  history  of  drag  coefficient  and  value  of  drag  near  convergent 

state. 
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Figure  6.5  Drag  coefficient  after  vortex  shedding  for  both  inlet  conditions. 
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Figure  6.6  Convergence  history  of  drag  and  lift  coefficient  using  bounce-back  inlet 

condition. 
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Figure  6.7  Convergence  history  of  drag  and  lift  coefficient  using  local  inlet  condition. 
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Figure  6.8  Convergence  history  of  Drag  coefficient  at  Re=2000  using  bounce-back 
condition  at  inlet,  upper  and  lower  boundaries. 
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Figure  6.9  Convergence  history  of  Drag  coefficient  at  Re=2000  using  local  inlet 
condition  at  inlet,  upper  and  lower  boundaries 
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Figure  6.10  Density  contour  at  time  equaling  to  30000  lattice  unit  at  Re=2000  using 
bounce-back  condition  at  inlet,  upper  and  lower  boundaries. 


Figure  6.1 1 Density  contour  at  time  equaling  to  30000  lattice  unit  at  Re=2000  using  local 
inlet  condition  at  inlet,  upper  and  lower  boundaries 
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Y 

Figure  6.12  Velocity  profile  along  y-  direction  at  nodes  which  are  0.2  lattice  unit  to  the 

right  of  inlet. 
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Figure  6.13  Shear  stress  profile  along  y-direction  at  nodes  which  are  0.2  lattice  unit  to  the 

right  of  inlet. 
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6.2  A Unified  Treatment  for  Solid  Wall  Boundary  Condition  in  LBE 

Recent  development  in  solid  wall  boundary  conditions  can  be  found  in  Filippova  and 
Hanel  (1998),  Mei  et  al.  (1999),  and  Bouzidi  et  al.  (2001).  In  order  to  establish  the 
necessary  background,  these  three  treatments  will  be  reviewed  in  the  following. 

The  lattice  nodes  on  the  solid  and  fluid  side  are  denoted  as  xb  and  Xf  respectively,  in 
Figure  6.14.  The  filled  small  circles  on  the  boundary,  xw,  denote  the  intersections  of  the 
wall  with  various  lattice  links.  The  boundary  velocity  at  xw  is  uw.  The  fraction  of  an 
intersected  link  in  the  fluid  region  is  A,  that  is: 


Obviously,  0<zl<l  and  the  horizontal  or  vertical  distance  between  xb  and  xw  is  Ack  on 
the  square  lattice.  Suppose  the  particle  momentum  moving  from  x/  to  xb  is  ea  and  the 
reversed  one  from  xb  to  X/  is  e u = -ea.  Denoting  the  post-collosion  properties  using  ~ , it 

is  seen  that  after  the  collision  step,  fa  on  the  fluid  side  is  known,  but  not  on  the  solid 
side.  Using  e ^ and  to  denote  the  velocity  and  the  distribution  function  coming  from 
a solid  node  to  a fluid  node,  it  is  noted  that  is  the  unknown  to  be  determined.  That  is, 

to  finish  the  subsequent  streaming  step,  fu{xb,t)  is  needed.  To  construct  f^(xb,t) 

based  upon  some  known  information  in  the  surrounding,  Filippova  and  Hanel  (1998) 
proposed  to  use  the  following  linear  interpolation: 


(6.15) 


fa  (xb’t)  -('\-Z)7a(xf’t)+Zfa\xb,  0 +2Wap 


3 


(6.16) 


with 


2c4 


9 


(e«-  uf)2  - 


y Uf  Uf] 

2 c2 


(6.17) 
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In  the  above,  Uf=  u(Xf,  t)  is  the  fluid  velocity  near  the  wall  and  «*/is  to  be  chosen. 


Ubf=  (A  -X)Uj! A + uJA  and 

Z = (2A-\)/t 

for 

A>\/2  (6.18a) 

Ubf-Uf  and 

Z = (2A-1)/(t-1) 

for 

A<\!2.  (6.18b) 

Mei  et  al.  (1999)  suggested  using  different  nodes  to  obtain  fa  when  zl<0.5  to  improve 
the  numerical  stability  over  the  scheme  of  Filippova  and  Hanel  for  the  channel  flow  when 
r is  near  0.5.  Bouzidi  et  al.  (2001)  presented  a simpler  boundary  condition  based  on 
bounce-back  at  arbitrary  position.  In  their  work,  a linear  scheme  and  a quadratic  scheme 
are  given  to  obtain  fa  at  node  inside  the  solid  wall.  The  linear  version  is  given  as 
following: 

+ 1) = 2A/0(r/  +e0)  + (l  -2A)/a(r/)  A<i  (6.19a) 

/*('/.'  + 0 = ^/>/+O  + ^/s(^+*.-)  A<i  (6.19b) 

They  also  reported  numerical  results  for  a circular  Couette  flow  and  flow  over  a 
periodical  array  of  cylinders. 

The  above  three  boundary  condition  treatments  all  have  second-order  accuracy  for 
curved  boundary.  The  difference  is  that  the  first  two  need  to  construct  a fictitious  fluid 
point  inside  the  solid  wall,  and  perform  a collision  step  at  that  node,  while  the  scheme  of 
Bouzidi  et  al  only  requires  the  known  values  of fa  on  the  fluid  side.  It  is  emphasized  that 
three  methods  need  to  treat  the  boundary  condition  separately  for  zl<0.5  and  zl>0.5. 

In  this  effort,  a unified  scheme  for  curved  boundary  condition  is  presented,  and  tested 
along  with  those  by  Filippova  and  Flanel  (1998),  and  Bouzidi  et  al.  (2000)  for  several 
fluid  flow  problems:  2-D  channel  flow  with  constant  pressure  gradient,  the  Stokes  first 
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problem,  and  uniform  flow  over  a column  of  circular  cylinder.  Finally  the  results  of  flow 
over  an  oscillation  plate  are  given. 

6.2.1  The  Unified  Boundary  Condition  for  Solid  Wall  Boundary 

When  simulating  the  evolution  of  fa  with  the  collision  and  streaming  steps,  one  needs 
to  address  several  issues.  For  example,  to  finish  the  streaming  step,/a  at  solid  boundary 
node  Xb  must  be  specified.  Usually  we  consider  obtaining  the  required  boundary  condition 
based  on  the  value  of  fa  after  the  collision  step.  We  can  also  perform  a streaming  step 
with  a null  boundary  value  fa,  and  then  add  the  null  boundary  value  with  appropriate  f„ 
after  the  streaming.  There  is  a difference  between  these  two  approaches,  because  the 
streaming  redistributes  fa  in  space.  Referring  to  Figure  6.14,  after  streaming,^  at  the 
fluid  node  f is  needed  while  the  values  of  f,s  and  fj$  are  known.  Using  a linear 
interpolation,  fw#  can  be  easily  found  as: 

/,8  = f ft  + (As  - f ft ) * A (6.20) 

Here  fw8  is  the  distribution  function  at  the  point  on  the  wall  where  the  link  along  e8 
direction  intersects  with  the  solid  wall.  To  ensure  the  no-slip  boundary  condition  (uw  = 0 
here)  on  the  wall,  considering  the  momentum  balance  in  each  direction,  we  set 

fw4  ~ fwi  (6.21) 

Using  the  /A  and  fj/4  (which  is  fjj4  at  the  fluid  node  ff ),ff4  can  be  obtained  using  a linear 
interpolation: 

ff  4 = /„ 4 + (Jjf  4 " L 4 ) /(l  + A)  * A (6.22) 

This  simple,  unified  formula  for  the  boundary  condition  is  valid  for  both  A > 0.5  and 


A <0.5. 
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The  extension  can  be  easily  made  using  a high  order  interpolation.  Take  the  second 
order  formulation  for  example. 

Let: 

Ti=y*8’  y 2 ~ f 

x,  = 1 ; x,  = 2 ; x3  =3  (6.23) 

x = 1 - A 


The  following  three-point  Lagrange  interpolation  is  used  to  obtain  fw$  : 


f.s  = £ 


k=\ 


x - X , 


J- 1 Xk  ~XJ 


n 

7=1 

\j*k 


■yk 


By  requiring 

f *4  f .i'8 

the  no-slip  boundary  condition  on  the  wall  is  ensured. 

Letting: 

y 1 = fw 4’  T2  =/^4’  T3  = /®r4 
x,=l  ; xi=2  + A;x3=3  + A 
x = A 


(6.24) 


(6.25) 


(6.26) 


and  employing  three  point  Lagrange  interpolation  again,  the  distribution  function  /)  at 
node  / can  be  readily  obtained.  In  Eq.  (6.26),  ‘‘^4  ” donates  the  node  which  is  one  lattice 
unit  away  from  the  node  ff  along  the  ^-direction. 

For  non-zero  wall  velocity,  uw,  it  can  be  easily  incorporated  into  the  present  unified 
formulation  by  adding  an  additional  momentum  to  fws.  given  in  (6.25) 

3 

La  =La+2WaP  — ea-U*’ 

C 


for  cc=  4 


(6.27) 
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where  p is  the  fluid  density  at  the  wall  which  can  be  obtained  by  using  the  extrapolation 
from  the  nearby  fluid  nodes.  To  avoid  computational  instability,  we  simply  s z\p  = p f for 

nearly  incompressible  flow. 

6.2.2  Computational  Assessment 

The  flow  over  a column  of  circular  cylinders  is  the  case  used  to  assess  the  impact  of 
the  boundary  treatment  on  the  accuracy  of  the  flow  field  around  a curved  boundary.  The 
flow  over  an  oscillating  zero-thickness  plate  is  used  to  assess  the  temporal  continuity  of 
the  solution  after  the  application  of  the  presently  proposed  boundary  treatment. 

6.2.2. 1 Pressure  driven  channel  flows 

The  grid  structure  for  the  2-D  channel  computation  is  shown  in  Figure  4.3.  A constant 
pressure  gradient  Vp  along  the  x-direction  is  applied  as  a body  force,  which  is  added  after 
the  collision  step  (He  et  al.  1 997): 


where  x is  the  unit  vector  along  the  x-direction.  At  steady  state,  the  exact  solution  for  the 
jc-velocity  profile  is 


where  H = Ny-3+2A  and  rj=y/H=(j- 2+  A)/H.  To  assess  the  computational  error  of  the 
LBE  solution  of  the  velocity,  u lbe  (y),  the  following  relative  Z-2-norm  error  is  defined 


(6.28) 


2 dx  pv 


(6.29) 


[Joh“Lc,0'MT'2 


exact 


exact 


(6.30) 
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In  the  LBE,  At=Ax=Ay=\.  Using  the  channel  height  H=Ny-3+2A,  the  dimensionless  grid 
size  (or  grid  resolution)  is  H~ 

The  periodic  boundary  condition  is  used  at  the  left  and  right  boundaries  of  the 
computational  domain.  At  the  lower  and  upper  walls,  three  wall  boundary  conditions  are 
used:  i)  Fillipova  and  Hand’s  method,  ii)  Bouzidi  et  a/’s  method  and  iii)  the  present 
method.  For  the  computations  reported  here,  Vp=  10'8,  r=0.51  with  double  precision. 

The  wall  slip  velocity  uw=ux(y=0)  is  evaluated  using  a second-order  extrapolation 
based  on  ux(y=A),  ux(y=\+A)  and  ux(y=2+A).  Flere  we  suppose  that  the  velocity  profile  of 
the  solution  is  parabolic  so  that  this  second  order  extrapolation  will  not  cause  additional 
extrapolation  error.  Since  the  true  wall  velocity  in  the  pressure  driven  channel  flow  is 
exactly  zero,  the  wall  slip  velocity  uw  provides  a measure  of  the  accuracy  for  the 
treatment  of  the  wall  velocity.  Figure  6.15  shows  the  dependence  of  wall  slip  velocity  to 
the  H (or  the  grid  resolution  FT1).  Here  uw  is  normalized  by  the  centerline  velocity 


Mmax 


Hz  dp 
8 pv  dx 


(6.31) 


The  second-order  convergence  of  uw  with  increasing  H is  observed  clearly  in  the 
Figure  6.15  for  zl=0.01,  0.5,  and  0.99.  It  is  worth  noting  that  although  both  the  present 
linear  and  quadratic  formulae  show  the  second-order  convergence,  the  magnitude  of  error 
of  the  quadratic  form  is  much  smaller. 

Figure  6.16  shows  the  dependence  of  the  relative  12-norm  error  on  the  channel  height 
H for  d=0.01,  0.5,  and  0.99.  The  second-order  accuracy  is  clearly  demonstrated  in  the 


range  of  H investigated.  It  has  been  well  established  that  the  accuracy  of  the  LBE 
method  for  the  interior  points  is  of  second  order.  The  fact  that  the  overall  accuracy  is  of 
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second  order  in  the  present  case  means  that  the  accuracy  in  the  boundary  condition  is  at 
least  of  second  order.  This  is  entirely  consistent  with  the  results  shown  in  Figure  6.15, 
which  shows  local  convergence  (y= 0). 

A careful  examination  of  the  relative  errors  in  Figure  6.15  and  Figure  6.16  indicates 
that  the  wall  slip  velocity  in  the  LBE  may  be  responsible  for  the  error  in  the  flow  field.  It 
is  thus  instructive  to  compare  the  shifted  velocity  profile, 

be  O')  = ui.be  O)  - (6.32) 

with  the  exact  velocity  profile  since  the  error  on  the  wall  in  the  shifted  profile  is  exactly 
zero  now.  The  relative  T2-norm  error  of  uLBE(y)  is 


E2  = 


{f  [ULBE(y)-^exacl(y)]2dy}U2 

[f«L  ,(ymU2 


(6.33) 


Figure  6.17  shows  the  result  of  Z,2-norm  error  using  the  Eq.  (6.33).  The  variation  of  E2  is 
non-monotonic  because  the  absolute  error  is  already  at  the  level  of  machine  accuracy.  In 
computation,  we  use  double  precision  and  uexac,  is  around  10‘4.  From  Figure  6.17,  it  is 
seen  that  LBE  does  not  generate  any  error  in  the  computation  for  the  channel  flow;  the 
error  is  entirely  from  the  boundary  treatment.  Hence  the  /.2-norm  error  given  by  the  Eq. 
(6.30)  can  be  used  to  characterize  the  error  associated  with  the  boundary  condition. 
Figure  6.18  shows  the  relative  error  by  Eq.  (6.31)  as  a function  of  Zl  using  the  present 
method,  FH’s  method,  and  method  of  Bouzidi  et  al.,  for  0<H<1.  It  can  be  seen  that  the 
present  linear  scheme,  linear  scheme  of  Bouzidi  et  al.,  and  FH’s  boundary  condition  have 
the  same  order  of  error  for  0<zl<l.  Bouzidi  et  al.' s quadratic  form  gives  a large  range  of 
error  for  0<H<1,  from  10'  to  10'  . The  present  quadratic  formula  yields  a more 
uniformly  distributed  error  in  the  whole  range  of  A.  It  is  interesting  to  note  that  the  error 
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goes  down  when  A increase  for  the  present  quadratic  formula,  which  is  opposite  to  the 
trends  shown  in  other  boundary  conditions.  The  present  linear  and  Bouzidi  linear 
boundary  conditions  have  very  similar  characteristics.  Unlike  the  FH  boundary 
condition,  the  present  boundary  conditions  suffers  no  computational  instability  for  r>  0.5 
and  any  value  of  A between  0 and  1.  Even  for  t=5.0+10~4,  the  present  boundary 
condition  gives  a correct  converged  result. 

6.2.22  Stokes  first  problem:  flow  due  to  an  impulsively  started  wall 

For  a wall  located  at  y=  0 that  is  impulsively  started,  an  unsteady  Stokes  layer  of 
thickness  0(4vt ) develops  near  the  wall.  For  a fixed-grid  computation,  the  error  at  small 
time  is  expected  to  be  large  due  to  insufficient  spatial  resolution.  In  the  LBE  method,  this 
is  also  compounded  by  the  use  of  fixed  St  (=Sx=Sy=\).  In  the  computation  the  wall 
velocity  is  set  as  V=0A.  Two  boundary  conditions  are  used  on  the  wall  in  the 
computation,  FFI  boundary  condition  and  present  boundary  condition.  The  relaxation 
time  z=0.7  gives  a kinematic  viscosity  v=0.067.  The  Stokes  layer  is  about  0.26  and  0.82 
< Sy=\  respectively  when  t=\.0  and  10  (lattice  unit).  Hence  there  is  not  enough  spatial 
resolution  when  t is  small  and  this  will  typically  cause  spatial  and  temporal  oscillation  in 
the  solution. 

Figure  6.19  shows  the  velocity  profiles  at  t=  9.5,  49.5,  and  99.5  (in  lattice  unit).  When 
t= 9.5,  the  velocity  has  large  relative  error  and  it  becomes  negative  when  y is  larger  than 
4.  The  oscillation  of  velocity  is  not  shown  in  the  Figure  6.19  because  its  absolute  value  is 
less  than  10'J.  Only  the  result  of  present  linear  boundary  treatment  is  shown  in  Figure 


6.19. 
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Figure  6.20a  and  Figure  6.20b  shows  the  temporal  variation  of  the  relative  Z,2-norm 
error  defined  as 


r _ (lo" lu LBE(y ,0  - « exact  CM)]2^}* /2  ^ ,,  ^ 

^2 — ~z T7~> (6.34) 

[%uLa(y.>m'2 

for  zl=0.01,  0.5,  and  0.99.  The  large  relative  errors  in  the  beginning  are  due  to  lack  of 
resolution  for  small  time.  It  should  be  emphasized  that  this  flow  at  small  time  is  difficult 
to  deal  with  for  any  computational  technique  due  to  the  singular  acceleration  and  large 
spatial  gradient.  It  is  interesting  to  note  that  the  present  quadratic  boundary  condition 
treatment  converge  to  the  same  magnitude  of  error  while  the  FH  and  the  present  linear 
boundary  formulas  diverge  at  a different  level.  In  such  a transient  flow,  the 
computational  accuracy  in  the  near-wall  region  is  typically  dictated  by  the  near-wall 
spatial  resolution  which  must  be  smaller  than  the  Stokes  layer  thickness  in  order  to 
resolve  the  local  flow  field.  In  a finite  difference  calculation  for  such  a flow,  St  and  Sx 
can  be  independently  chosen.  If  Sx  is  not  sufficiently  small,  further  reduction  in  St  will 
not  lead  to  improvement  in  accuracy.  As  the  Stokes  layer  grows  to  a certain  thickness, 
the  spatial  resolution  becomes  adequate  and  the  accuracy  then  improves. 

6.2.2.3  Steady  uniform  flow  over  a column  of  cylinders 

For  a uniform  flow  over  a column  of  circular  cylinders  of  radius  r with  center-to- 
center  distance  denoted  by  //,  symmetry  conditions  for /a’s  are  imposed  at  y=  +H/2.  At 
the  inlet,  the  uniform  velocity,  u=V,  is  specified  using  FH  boundary  condition.  At  the 
exit,  a simple  extrapolation  is  used, 


MNx,j)  = 2MNx-\J)-MNx-2,j)  for  a=4,  5,  and  6. 


(6.35) 
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On  the  surface  of  the  circular  cylinder,  the  FH’s  and  the  present  methods  are  both 
used  to  update  the  fd s.  To  reduce  the  effect  of  the  inlet  disturbance,  the  cylinder  is 
placed  20  radii  downstream  of  the  inlet.  The  outlet  is  placed  35  radii  to  the  right  of 
center  of  cylinder.  In  the  computation  z=0.52,  w,„/e,=0.0513.  Re=100  based  on  the 

diameter  of  the  cylinder  and  inlet  velocity.  The  drag  coefficient  is  calculated  using  the 
momentum  exchange  method  for  both  boundary  condition  at  two  different  cylinder  center 
positions  (130.0,  65.0)  and  (130.2,  65.0).  Different  positions  of  the  cylinder  center  give 
different  A configurations  in  the  boundary  description. 

Table  6.1  shows  the  results  of  drag  coefficient  using  the  FH,  the  present  linear,  and 
the  present  quadratic  formulae.  The  present  quadratic  formula  gives  a closer  value  in 
comparison  with  the  value  of  1.248  given  by  Fornberg  (1991)  for  both  cylinder  positions. 
Figure  6.21  shows  the  centerline  velocity  variations  before  and  after  the  cylinder  obtained 
using  two  kinds  of  boundary  conditions  for  the  case  with  the  center  of  the  cylinder  at 
(130.0,65).  The  present  linear  form  and  FH’s  method  give  very  similar  velocity  profiles 
in  regions  with  sharp  gradient  near  the  front  stagnation  point,  very  similar  values  for  the 
length  of  the  separation  bubble,  the  maximum  of  the  separation  bubble  velocity,  and  the 
recovery  of  the  wake  velocity.  This  result  is  expected  because  the  results  of  Co  of  the  two 
cases  are  very  close.  The  present  quadratic  formula  gives  slightly  different  length  of  the 
separation  bubble,  and  hence  the  recovery  of  the  wake  velocity. 


Table  6.1  Comparison  of  Drag  coefficient  using  three  kinds  boundary 


X 

Present  Linear 

Present  quadratic 

FH 

CD, 

(Fornberg,  1991) 

130.0 

1.275 

1.247 

1.271 

1.248 

130.2 

1.290 

1.251 

1.274 
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6.2.2.4  Flow  over  a vertically  oscillating  zero-thickness  flat  plate 

The  computational  domain  is  shown  in  Figure  6.22.  The  inlet,  lower,  upper,  and  outlet 
boundary  condition  are  set  to  be  the  same  as  those  for  flow  over  a cylinder  discusses  in 
the  above  section.  There  are  40  grids  along  the  plate.  To  avoid  the  crossing  of  the  nodes 
by  the  plate,  the  plate  is  moving  according  to: 

y Plate  = y middle  + 0.49  • sin(f  / 1 000) 

In  the  computation,  7=0.52;  Re= 20  based  on  the  length  of  plate;  =0.0083.  FH’s 
method.  Bouzidi  et  aV s form,  and  present  linear  from  boundary  conditions  are  used  on 
the  surface  of  the  plate.  Figure  6.23a  and  Figure  6.23b  show  the  variation  of  drag  when 
plate  moves  up  and  down.  We  can  see  that  there  is  a discontinuity  in  the  results  using 
FFFs  and  Bouzidi  et  aV s formulae  when  the  plate  crosses  the  middle  of  the  lattice.  The 
present  boundary  condition  gives  a treatment  results  in  a smooth  temporal  variation. 

6.2.3  Conclusion 

In  this  work  a second-order  accurate  boundary  condition  treatment  for  the  lattice 
Boltzmann  equation  is  proposed.  A series  of  studies  are  conducted  to  systematically 
validate  the  accuracy  and  examine  the  robustness  of  the  proposed  boundary  condition  in 
steady  and  unsteady  flows  involving  flat  and  curved  walls.  Compared  with  the  existing 
method  for  treating  boundary  condition  in  the  lattice  Boltzmann  method,  the  proposed 
treatment  has  the  following  advantages,  (i)  The  scheme  is  unified  for  A>0.5  and  A<0.5 
and  no  discontinuity  in  boundary  treatment,  (ii)  The  boundary  treatment  is  very  robust 
and  surfer  no  computational  instabilities  in  the  pressure  driven  channel  flow  for  r>  0.5. 
(iii)  Compared  with  the  other  curved  boundary  condition,  the  present  quadratic  condition 
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gives  overall  better  result  under  otherwise  identical  computational  parameters,  (iv)  The 
present  scheme  is  very  simple  and  easily  extendable. 


fff@ 

\ 


Figure  6.14  Layout  of  the  regularly  spaced  lattices  and  curved  wall  boundary. 
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Figure  6.15  Quadratic  convergence  of  the  wall  slip  velocity  for  present  boundary 
condition  in  constant  pressure  driven  channel  flow  . 


A 


Figure  6.16  Dependence  of  relative  Z,2-norm  error  using  Eq.  (6.30)  for  the  present 
boundary  conditions  on  the  lattice  resolution  H in  steady  state  pressure-driven 
channel  flow  simulations. 


Relative  L,  norm  error 
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H 

Figure  6.17  Dependence  of  relative  T^-norm  error  using  Eq.  (6.33)  for  the  present 
boundary  conditions  on  the  lattice  resolution  H. 


Figure  6.18  Relative  Z,2-norm  error  using  Eq.  (6.30)  as  a function  of  A for 
different  boundary  conditions  in  steady  state  pressure-driven  channel  flow 
simulations. 
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Figure  6.19  Velocity  profiles  at  different  time  in  an  impulsively  started  plate  (A= 
Vi)  using  present  linear  boundary  condition. 


Figure  6.20a  Relative  Z.2-norm  error  of  the  velocity  profile  ux(y)  during  the 
initial  transient  of  the  impulsively  started  plate  with  various  values  of  A for  the 
present  boundary  conditions. 
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Figure  6.20b  Relative  Z,2-norm  error  of  the  velocity  profile  ux{y)  during  the 
initial  transient  of  the  impulsively  started  plate  with  various  values  of  A for  the 
present  quadratic  and  FH  boundary  conditions. 


Figure  6.21  Centerline  velocity  variation  for  a uniform  flow  over  a column  of 
cylinders  with  FH,  present  quadratic,  and  present  linear  boundary  conditions. 
Center  of  cylinder  is  at  (130,65)  lattice  unit. 
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Figure  6.22  Computational  domain  of  for  flow  over  oscillating  zero-thickness  flat 
plate. 


FH' s method 


Figure  6.23a  Drag  coefficient  Cd  for  flow  over  oscillating  zero-thickness  flat  plate 
using  three  different  boundary  conditions:  FH.  Bouzidi  et  al.  linear,  and  present 
linear  scheme. 
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Figure  6.23b  Enlarged  view  of  Cd  in  the  circled  region  in  Figure  6.22a  when  plate 
crosses  the  middle  of  grids.  The  figure  demonstrates  the  discontinuity  in  the 
temporal  variation  caused  by  the  FH  and  Bouzidi  et  al ’s  boundary  treatment. 
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6.3  An  Improved  Moving  Wall  Boundary  Condition  in  LBE 
6.3.1  Issues  in  the  Moving  Boundary  Problems 

Unlike  the  problems  for  flow  over  a stationary  object,  where  the  position  and 
configuration  of  the  boundary  do  not  change  with  time,  moving  boundary  problems 
generate  additional  difficulties  during  the  course  of  the  computation.  To  facilitate  the 
discussion,  a number  of  terminologies  are  introduced  first.  Here  we  first  give  some 
definitions  for  better  description.  A fluid  boundary  node  (FBN)  is  a node  neighboring  the 
wall  and  residing  the  fluid.  A solid  boundary  node  (SBN)  is  the  node  neighboring  the 
wall  and  inside  the  solid.  A pre-fluid  boundary  node  (PFBN)  is  a fluid  node  neighboring 
a fluid  boundary  node  in  one  direction  and  will  become  a fluid  boundary  node  in  the 
future;  see  Figure  6.24  for  details.  When  a boundary  moves,  there  are  two  possible 
changes  in  the  grid  arrangement.  The  first  is  that  a PFBN  becomes  a new  FBN  and  a 
FBN  becomes  a new  SBN.  The  second  is  that  a SBN  becomes  a new  FBN.  In  the 
second  case  information  at  the  newly  created  FBN  is  not  available  directly.  Suppose  that 
wall  velocity  is  small,  when  a PFBN  becomes  a FBN  and  directly  exchange  momentum 
with  the  wall,  the  physical  variables  at  that  new  FBN  node  should  not  have  a sudden 
jump. 

To  illustrate  the  issue  related  to  a possible  sudden  jump  in  flow  variables,  a one- 
dimensional illustration  is  shown  in  Figure  6.25.  At  time  t,  the  boundary  is  at  position  1. 
After  one  time  step  (t+1),  the  boundary  is  at  position  2.  In  this  evolution,  the  PFBN  A 
becomes  a FBN  and  the  FBN  B becomes  a SBN.  Supposed  that  the  wall  velocity  is 
small,  the  distance  from  the  PFBN  A to  the  solid  wall  at  time  t and  the  distance  from 
FBN  A to  the  solid  wall  at  time  t+1  should  only  have  a very  small  change,  so  the 
physical  variables  at  the  PFBN  A and  FBN  A must  change  smoothly.  That  means  any 
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sudden  jumps  of  the  variables  at  the  node  A during  this  procedure  are  nonphysical.  This 
is  the  first  problem  to  be  addressed  in  the  moving  boundary  problem.  If  the  bounce-back 
condition  is  used,  one  would  not  be  surprised  to  see  that  the  physical  variables  at  node  A 
exhibits  a discontinuity  when  the  solid  boundary  crosses  a nodal  location.  The  bounce- 
back  scheme  handles  the  curved  wall  as  zig-zag  steps.  When  the  boundary  crosses  the 
node  B,  the  real  position  of  boundary  will  move  from  the  middle  of  the  nodes  B and  C to 
the  middle  of  the  nodes  A and  B.  This  sudden  jump  of  boundary  position  causes  a large 
jump  of  variables  at  the  node  A. 

This  kind  of  jump  is  not  expected  to  arise  when  the  integrity  of  curved  boundary  is 
honored  to  the  second  order  accuracy  for  0<A<7.  In  what  follows  the  FH’s  boundary 
condition  is  analyzed  to  further  illustrate  how  large  the  change  in  physical  variables  at 
node  A would  results  from  this  procedure.  For  smooth  transition,  it  is  needed  to 

maintain /5(xfl)?+y  « /5  (■*«),  which  will  be  streamed  into  node  A at  t+1  and  t 
respectively. 

At  time  t,  the  node  B is  the  FBN  and  A « 0.  Using  FH’s  boundary  condition 
treatment,  see  Eq.  (6.16): 

~ 3 

f 5{xc ,t - \)=(\-x)  f\{x B,t  - \)  + xf\  {xBit-\)+2w5Pw—  e5  uw  (6.36) 

c 

because  of  A « 0,  Using  Eqs.  (6.17)  and  (6.18)one  has 

= (6.37) 

X = ~ l/(r-l)  (6.38) 

Substituting  Eq.  (6.37)  and  (6.38)  into  Eq.  (6.36),  one  obtains 
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7 (*c  >t  - 1)  = /,  (xA  ,t  - 2)  +2w5pw  e5  • «„ 


(6.39) 


The  steaming  step  gives 


f5(xB,t)  = f5(xc,t-\) 


(6.40) 


After  collision,  f5(xB,t)  becomes: 


r r 


Substituting  Eq.  (6.40)  into  Eq.  (6.41),  one  obtains: 


(6.41) 


(6.42) 


7(*«>o=  75(Acc,t-2)-(i-7)+//?(xfl,o-- 

r r 


The  equilibrium  distribution  function  in  Eq.  (6.42)  can  be  determined  as 

faq(xg,t)=wap(xB,  0 [ l + —ea-  ub  (0+~r(ea  uB{t)f  —2—  (, uB  (t))2  (6.43) 

c 2c  2c 

where  a=5.  Using  equations  (6.36)-(6.43)  and  the  condition  that  uB{i)~uw,  since  A ~ 0 
one  obtains: 


fs(xH,t)  =(f\(xA,t -2)+2wspw  —e5-uw)-(  1 — ) + 


T 

T 


3 9 2 3 

w5p(xB,  t)  [1  + —e$  uw  + — -(ef  «,„)  — uw  ■ uw] 


c2  2 c4  2 c2  t 


(6.44) 


At  time  t+1,  the  node  A is  the  new  FBN,  and  A changes  from  zl«0  to  A&l.  From 
FFTs  boundary  condition  one  obtains: 

f5(x/j,t)  = fl(xA,t)-(l--)  + fr)(xH,t)--+2wspw  uw  (6.45) 


and 
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fP(xB,t)=w,p(xA,  t)  [1  + \er  Uw  +-\(c7-  uA( t))2 7(uA  (t))2  ] (6.46) 

c2  2c4  2c2 

Using  equations  (6.45)  and  (6.46)  one  obtains: 

1 3 

fs(xH,t)  = fl(xA,t)-(l--)+(w,p(xA,  t)  [1  + —er  uw 

T C1 

uA(t))2  — — -{uA  (t))2])  --+2w5p  ~ e5  ■ uw  (6.47) 
2c4  2c 2 t c2 

Subtracting  Eq.  (6.44)  from  (6.47),  one  obtains: 

f 5 (XB  + 1)  ~ (Xft  ; 0 — [,/|  (XA  ’ 0 — f\  (XA  T ~ 2)](1  ) $p w — c5  • Uw+ 

T C 

3 9 2 3 

Wp\Xiu  f)  [1  + — C s'  Uw  t"  —{pf  U\: ) — Uw  ' U,v\- 

c2  2c4  2c2 

w,p(xA,  t+1 ) [1  + \er  Uw  + —^-(c/  UA( t+1))2  — ~(uA  (t+1))2]}  •-  (6.48) 

c2  2c4  2c2  t 

In  the  ideal  case,  as  the  solid  moves  crossing  the  node  B,  variables  at  position  A should 

not  experience  a sudden  jump.  This  implies  that  the  LHS  terms  are  small  and  would 

require  the  following  conditions  in  Eq.  (6.48): 

pixA,t+l)&  p(xB,t)  ; uA(t+l)*uw  (6.49) 

In  the  actual  computation,  these  conditions  cannot  be  satisfied  because  variables 

appearing  in  Eq.  (6.49)  belong  to  the  different  grids  at  time  t and  t+1. 

The  second  problem  in  moving  boundary  problems  is  about  the  creation  of  new  fluid 

nodes  (see  Figure  6.26).  The  SBN  at  the  location  E at  time  t becomes  a new  FBN  at  time 

t+1.  At  time  t+1,  the  information  for  f’s  at  the  node  E does  not  exist  and  need  to  be 

specified  using  the  information  from  the  nodes  around  the  new  FBN  E.  There  are 

different  ways  to  specify  the  distribution  function  at  the  node  E.  Here  we  first  use  a 
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linear  extrapolation  to  give  all  distribution  functions  at  node  D.  An  additional  step  is  then 
taken  to  reduce  the  possible  oscillation  caused  by  this  linear  extrapolation.  This 
additional  step  can  be  described  as  placing  an  additional  collision  in  the  fluid  nodes 

neighboring  the  new  FBN  and  then  streaming  fa  from  the  neighboring  nodes  to  the  new 
FBN  E.  Then  a boundary  condition  is  needed  to  supply  , where  eu  = -ea . In  the 
second  step  only  the  value  of  fa  at  the  new  FBN  is  changed.  Although  a collision  step 
has  been  added  at  the  nodes  neighboring  the  new  FBN,  the  values  of  fa  at  those  nodes  do 

not  change.  The  post-collision  values  are  used  only  to  stream  to  the  new  FBN  E. 

The  third  problem  appearing  in  the  moving  boundary  problem  is  the  way  the 
singularity  is  handled.  For  the  stationary  boundary  condition,  the  ability  to  handle  the 
spatial  and  temporal  accuracy  and  geometric  singularity  is  very  important.  With  the 
geometric  singularity  such  as  a sharp  convex  corners,  the  solution  may  posses  spatial 
oscillation  and  cause  computational  instability  as  well  as  reduced  accuracy  of  the 
solution.  In  the  moving  boundary  problem,  the  geometric  singularity  causes  even  a 
bigger  problem  because  the  moving  boundary  changes  the  position  of  the  singularity  and 
redistributes  the  A in  space. 

6.3.2  Formulation  for  the  Moving  Boundary  Condition 

In  LBE,  the  evolution  of  fa  ’s  involves  collision  and  streaming  steps.  To  complete 

the  streaming  step,  fa  at  solid  boundary  node  x*  in  Figure  6.24  must  be  specified. 

Referring  to  Figure  6.24,  after  the  streaming  step,  f/4  at  the  fluid  node  f needs  to  be 
assigned  while  fb&  andjfc  are  known.  Using  linear  interpolation,  fw8  can  be  found  as: 


fwS  ~ f f%  + (As  f f% ) * A 


(6.50) 
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where  fws  is  the  distribution  function  at  the  point  on  the  wall  and  the  link  along  the  eg 
direction  intersects  with  the  solid  wall.  To  ensure  the  no-slip  boundary  condition  at  w,  an 
additional  term  is  need  to  account  for  the  momentum  exchange  with  the  solid: 

/„4  =/ws  +2w4pw^4  uw  (6.51) 

c 

where  /},.  is  the  fluid  density  at  the  wall,  which  can  be  obtained  using  a suitable 
extrapolation. 

The  key  step  in  designing  this  moving  boundary  condition  is  to  add  an  additional 
collision  step  corresponding  to  the  boundary  velocity  «w  on  the  wall  to  obtain  /w4  and 
then  to  stream  it  by  one  lattice  in  the  direction  of  e*. 

7. 4 =/„. -a--) +/;?•-  (6-52> 

T T 

where 

/* 4 = [1  + wH+  — — (ej  «„,)2  — %r  (uw  )2]  (6.53) 

c2  2c4  2c2 

Based  on  the  assumption  of  the  slow  flow,  an  interpolation  is  used  to  obtain  to  //4  as, 

//.=/,, +(7,. -/.<)•  A (6.54) 

In  this  approach,  the  boundary  condition  does  not  need  to  be  treated  differently  for 
A > 0.5  and  A < 0.5 . However  this  additional  collision  step  is  the  key  difference 
between  this  scheme  and  the  one  described  in  Section  6.2. 

6.3.3  An  Analysis  on  the  Pressure  Change  near  the  Wall  and  Momentum 
Exchange  between  Fluid  and  Solid 

When  a boundary  crosses  a nodal  location,  it  is  desirable  to  keep  the  physical  variables 
at  the  new  FBN  and  the  momentum  flux  between  the  fluid  and  the  solid  continuous.  To 
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illustrate  that  the  present  boundary  condition  can  meet  both  requirements,  the  one 
dimension  example  is  again  used. 

First  we  show  that  physical  variables  at  new  FBN  can  be  maintained  continuous  (see 
Figure  6.25)  when  the  solid  boundary  crosses  a node.  As  shown  in  Section  6.3.1,  f5(xB) 
needs  to  be  examined  to  see  if  the  temporal  changes  of  variables  at  node  A are  smooth. 
At  time  t,  the  node  B is  the  FBN  and  A»0.  Using  t he  present  boundary  condition,  from 
Eq.  (6.50)  one  has: 

MxB,t)  (6.55) 

From  Eq.  (6.51),  one  has 

3 

f5(xw,t)  = f](xB,t)  + 2wap—ea  uw  (6.56) 

c 

Equation  (6.54)  gives 

3 

J5(xH^)  = MxA,t-\)  + 2wap~ea  uw  (6.57) 

c 

After  one  collision  step: 

75(**,0  =fs(xB,t)  -(1-7)  + fseq(xH,t)-T 

T 

= (MXA’t~'[)  + 2waP-Tea  •O-(l--)  + //9(x/j,0-r  (6'5g) 

c r 

The  equilibrium  distribution  in  Eq.  (6.58)  is  determined  as: 

fa‘'(xB,t)=wap(xB,t)[\+^Tea  ■uB(t)  + ^T(ea  uB(t))2  -^(«„(0)2]  (6.59) 

c 2c  2c 

At  time  t+1,  the  node  A is  the  new  FBN  and  1.  Using  the  present  boundary  condition, 


one  obtains: 
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f5(xw,t  + \)  = fi(xA,t  + l)  + 2wapw-ea-uw  (6.60) 

c 

After  one  collision  step  at  the  wall,  we  obtain: 

fs(xw,t+ 1)  =f5(xlv,t  + l)-(\--)  + f5ev(xw,t  + \)-T 

T 

=fi(xA,t)  + 2wap~ea-uw  -(1  --)+f5eq(xB,t  + \)-r  (6.61) 

C T 

where  the  equilibrium  distribution  function  is  given  by: 

3 9 3 

f?(xw,t  + 1)=  wap(xw,t  + 1)[1+  — ea  -uw  +— r(eB  -uw)2  - — yuw  -uw]  (6.62) 

c 2c  2c 

Subtracting  Eq.  (6.58)  from  (6.61),  one  obtains 

75(JCJJ,/+1)-75(XJ,,/)=1^  (x^  , r + 1) + yj  (x^ , / — i)  J • (l  - — ) 

T 

+ [Aeq(xw,t  + \)-f5“’(xB,t)\-T  (6.63) 

Since  xB  » jc w,  uB(t)&uw  and  p(xw,t + \)&p(xB,t) , one  has 

f5(xB,t  + l)*f5(xB,t)  (6.64) 

This  shows  that  the  present  boundary  condition  can  ensure  smooth  variation  of  f’s  at  node 
A. 

In  LBE,  the  momentum  exchange  is  used  to  find  the  force  acted  on  the  object.  The 
momentum  exchange  should  not  experience  an  artificially  large  fluctuation  when  a solid 
boundary  crosses  anodal  line.  Indeed,  as  will  be  presented  next,  the  present  boundary 
condition  can  satisfy  this  requirement  as  well.  In  the  derivation  the  node  A is  used  as  the 
reference  point  (see  Figure  6.25).  At  time  t,  the  momentum  exchange  Fs(t)  in  directions 


e\  and  e$  is: 
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F5(t)  = fl(xB,t)  + fx(xA,t-\)+  2w5p\e5-uw  (6.65) 

c 

where  /, (xB,t) is  determined  by: 

7i  (*B  ,t)  = fi(xA,t-l)-(l~-)  + fiq  (xB,t)  • t (6.66) 

r 

Substituting  Eq.  (6.66)  into  (6.65),  we  obtain: 

F5(t)  = fx(xA,t-\)  + fx(xA,t-\)-(\--)  + fxq(xB,t)-r  + 2w5p\e5 -uw  (6.67) 

r c 

At  time  t+1,  the  momentum  exchange  is: 

F5(t  + 1)  = f(xA,t  + 1)  + 

(f\(xA,t  + \)  + 2 w5p\e5  ■ uw)-(\-—)  + fsq (xw,t  + 1) • r (6.68) 

C T 

Adding  fxq(xw,t  + 1)-  z on  both  sides  of  Eq.  (6.68),  we  obtain 
Fs(t  + 1)  = f](xA,t  + 1)  + 

fx(xA,t  + 1)  ■(\--)  + fxq(xw,t  + \)-T  + 2w5p\e5  -uw  (6.69) 

T C 

Subtracting  Eq.  (6.67)  from  (6.69),  one  obtains, 

F5(t  + l)-F5(t)=[fx(xA,t  + \)-f](xA,t-\)\2--) 

T 

+fr(x„,t+\)-T-fr(xB,t)-T  (6.70) 

Since  feq (xw  J + 1)*  fxeq(xB,t)  and  fx(xA,t  + \)&fx(xA,t-\),  one  obtains 

F5 (/  + 1) ~ F5(t)  (6.71) 

This  demonstrates  that  the  momentum  exchange  can  be  maintained  continuous  from  t to 
t+1  when  a lattice  node  is  crossed  by  the  solid  boundary. 
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6.3.4  Computational  Assessment 

In  what  follows,  we  consider  flow  problems  without  singularity  to  test  the  accuracy  of 
the  present  boundary  condition.  First  we  consider  a channel  flow  under  constant  pressure 
gradients  with  analytic  solutions  to  assess  the  spatial  accuracy.  In  this  case  the  wall 
velocity  is  zero.  Flow  over  impulsively  started  cylinder  and  flow  over  stationary  cylinder 
in  a channel  are  then  used  to  test  the  accuracy  of  moving  boundary  treatments  and 
Galilean  invariance.  Finally  flow  due  to  a vertically  oriented,  infinitely  long  plate 
moving  sinusoidally  in  the  x-direction  is  considered.  The  present  boundary  condition  and 
FH’s  boundary  condition  are  used  for  comparison. 

6.3.4. 1 Pressure  driven  channel  flows 

The  configuration  and  parameters  of  this  computation  are  the  same  as  those  of  the  case  in 
Section  6.2. 3.1.  The  second  order  convergence  of  uw  is  observed  for  zl=0.01,  0.5,  and 
0.99  in  Figure  6.27. 

6.3. 4.2  Flow  over  a impulsively  started  cylinder  in  a channel 

The  channel  has  the  dimension  Z,xJF=401xl41.  The  initial  position  of  the  cylinder 
center  is  at  the  middle  point  of  the  domain  in  the  y direction.  The  radius  of  cylinder  is 
15.5  lattice  unit.  At  t= 0,  the  cylinder  is  at  rest.  Then  the  cylinder  is  impulsively  started 
with  a constant  velocity  U= 0.00215  moving  in  the  negative  x-direction.  On  the  left, 
upper  and  lower  boundaries  the  velocity  is  set  to  be  zero.  On  the  right  side  boundary,  the 
zeroth  order  extrapolation  is  applied  for  fa  ’s.  With  t=0.52,  the  Reynolds  number  based 

on  the  diameter  of  the  cylinder  is  10.  When  the  distance  crossed  by  moving  cylinder 
exceeds  1 unit  of  lattice,  the  entire  grid  and  cylinder  is  moved  back  to  the  right  by  one 
unit  to  ensure  that  the  distance  from  the  left  side 
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remains  constant  with  a lattice.  This  configuration  of  the  flow  is  nearly  equivalent  to  the 
flow  with  the  following  initial  and  boundary  conditions:  the  fluid  and  cylinder  is  at  rest  at 
t=0;  then  the  fluid,  the  upper  wall  and  lower  wall  start  to  move  impulsively  with  velocity 
£7=0.00215  with  the  velocity  of  the  left  boundary  is  also  set  to  £7=0.00215.  However 
there  is  still  a major  computational  difference.  For  the  stationary  cylinder,  the  number  of 
solid  boundary  nodes  is  fixed.  For  the  moving  cylinder,  the  number  changes  with  time. 
Thus,  it  is  expected  that  there  will  be  larger  transient  variation  in  the  moving  cylinder 
case. 

Figure  6.28a  shows  the  history  of  the  drag  on  the  cylinder  for  both  cases.  It  is  seen 
that  the  transient  drag  with  two  boundary  treatments  agree  with  each  other  very  well. 
Figure  6.28b  shows  the  final  stage  of  computation  when  steady  state  is  expected.  The 
oscillation  caused  by  the  moving  boundary  is  about  4%  of  the  mean  value.  In  the  force 
evaluation,  the  momentum  exchange  is  used.  In  this  case,  it  is  estimated  that  the  method 
of  force  evaluation  and  the  introduction  of  the  new  fluid  nodes  are  responsible  for  most 
part  of  the  oscillation.  Figure  6.29a  and  Figure  6.29b  show  the  density  contour  in  the 
region  near  the  cylinder  with  the  density  being  normalized  by  the  largest  density  in  front 
of  the  cylinder.  The  same  density  level  is  set  for  both  cases.  The  oscillation  in  density 
contour  is  very  small  in  both  cases,  but  the  spatial  distributions  of  density  in  space  are  not 
exactly  the  same.  This  difference  is  attributed  to  the  change  of  the  number  of  the 
boundary  nodes  in  the  moving  cylinder  case. 

6.3.4.3  A flow  due  to  a sinusoidally  oscillating  Plate 

The  computational  domain  is  shown  in  Figure  6.30.  The  velocity  V=0.0  is  set  as  the 
initial  condition.  A periodic  boundary  condition  is  imposed  at  the  upper  and  lower 
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boundaries.  The  zeroth  order  extrapolation  is  used  at  the  left  and  right  boundaries.  The 
plate  is  moving  at  the  velocity: 

u=U*sin(6.28/\000*t),  £7=0.01  (6.72) 

The  height  of  plate  is  11  lattice  unit.  On  the  surface  of  plate,  the  FH’s  boundary 
condition  and  present  moving  wall  boundary  condition  are  used.  The  number  of  new 
FBNs  generated  by  moving  plate  equals  to  the  height  of  the  plate  in  the  lattice  unit. 

Figure  6.31  and  Figure  6.32  show  the  history  of  pressure  at  a point  in  the  fluid  field 
that  is  one  grid  away  to  the  left  of  the  plate  in  x-direction  using  two  different  boundary 
condition  treatments.  The  pressure  at  this  specified  point  is  obtained  by  linearly 
interpolating  the  values  from  the  FBN  and  PFBN  on  the  left  side  of  the  plate.  The 
present  boundary  condition  produces  little  artificial  small-scale  fluctuations  as  the  plate 
crosses  a nodal  line.  However,  FH’s  treatment  results  noticeable  small  scale,  large- 
amplitude  disturbances  in  the  temporal  pressure  variation.  It  is  interesting  to  note  that 
present  boundary  condition  can  also  reduce  the  pressure  jump  when  new  fluid  nodes 
generate. 

6.3.5  Conclusions 

Various  issues  associated  with  the  moving  boundary  problem  have  been  examined.  A 
new  boundary  condition,  which  can  ensure  the  continuity  of  the  flow  variable  as  well  as 
momentum  exchange  as  boundary  crosses  the  lattice  nodes,  is  presented.  Various  cases 
show  that  the  proposed  boundary  condition  can  greatly  reduce  the  oscillation  caused  by 
the  moving  boundary.  The  new  boundary  condition  is  easy  to  implement. 
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Figure  6.24  Layout  of  the  regularly  spaced  lattices  and  curved  wall  boundary  and 
definition  of  nodes. 
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Figure  6.25  The  generation  of  a new  FBN  and  SBN  in  a moving  boundary 
problem. 
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Figure  6.26  The  generation  of  a new  FBN  in  a moving  boundary  problem. 
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Figure  6.27  Quadratic  convergence  of  the  wall  slip  velocity  for  present  boundary 
condition  in  a constant  pressure  driven  channel  flow. 
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Figure  6.28a  Convergent  history  of  drag  coefficient. 
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Figure  6.28b  Drag  coefficient  at  the  final  stage  of  computation. 
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Figure  6.29a  Density  contour  for  the  flow  over  a stationary  cylinder. 
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Figure  6.29b  Density  contour  for  the  flow  over  a moving  cylinder. 


20 


10 


>- 


0 


U=0.01 

u=U*sin(6.28/1 000*t) 


(-)  (+) 


-10 


_l_ 

930 


_L_ 

940 


I ....  I ....  I ....  I 

950  960  970  980 

X 


Figure  6.30  Computational  domain  for  flow  due  to  a sinusoidally  oscillation  plate. 
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FH's  boundary  condition 


Figure  6.31a  History  of  pressure  at  a point  which  is  one  grid  away  to  the  left  of 
the  plate.  The  details  of  circle  regions  A and  B are  given  in  Figure  6.31b  and 
Figure  6.31c. 


FH's  boundary  condition 


Figure  6.31b  Enlarged  region  A in  Figure  6.31a 
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FITs  boundary  condition 


Figure  6.31c  Enlarged  region  B in  Figure  6.31a 
Figure  6.31  Pressure  characteristic  at  a point  which  is  one  grid  away  to  the  left  of 
the  plate.  FFFs  boundary  condition  is  used  in  computation. 


present  boundary  condition 


Figure  6.32a  History  of  pressure  at  a point  which  is  one  grid  away  to  the  left 
of  the  plate.  The  details  of  circle  regions  A and  B are  given  in  Figure  6.32b 
and  Figure  6.32c. 
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present  boundary  condition 


Figure  6.32b  Enlarged  region  A in  Figure  6.32a 


present  boundary  condition 


Figure  6.32c  Enlarged  region  B in  Figure  6.32a 

Figure  6.32  Pressure  characteristic  at  a point  which  is  one  grid  away  to  the  left  of 
the  plate.  Present  boundary  condition  is  used  in  computation. 


CHAPTER  7 

SUMMARY  AND  FUTURE  WORKS 

Compared  with  NS  equations,  the  LBE  method  has  its  own  advantages  and 
disadvantages.  The  disadvantages  of  LBE  method  include:  firstly  it  require  large  a mount 
of  memory  to  store  the  distribution  function;  secondly  solution  is  generally  dependent  on 
the  time  which  means  capturing  a steady  state  solution  need  a long  time.  The  advantages 
of  LBE  include:  firstly  scheme  is  simple;  secondly  evolution  rule  is  local  and 
computation  is  easy  to  be  parallelized. 

The  issues  needed  to  be  addressed  in  the  LBE  method  include: 

1 . Increasing  the  efficiency  of  computation 

2.  Implementing  LBE  in  high  Re  number  flows 

3.  Evaluating  the  force  on  a solid  body  in  the  LBE  simulation 

4.  Improving  the  boundary  treatments 

5.  Implementing  LBE  in  moving  boundary  problems 

6.  Implementing  LBE  in  turbulent  flows 

7.  implementing  LBE  in  compressible  flows 

In  this  dissertation,  most  effort  has  been  put  to  address  the  issues  in  the  first  five 
topics  mentioned  above.  This  effort  includes: 

1 . A multi-block  method  in  the  LBE  method 

To  increasing  the  computational  efficiency  in  the  LBE  method,  a multi-block 
method  is  developed.  In  this  multi-block  method: 
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• The  formulae  for  information  exchange  at  the  interface  between  coarse 
and  fine  blocks  are  derived  to  ensure  the  mass  and  momentum 
conservation,  and  stress  continuity  across  the  interface. 

• A customized  time  matching  between  coarse  and  fine  blocks  is  designed 
to  ensure  the  smooth  transferring  of  physical  variables  between  coarse  and 
fine  blocks. 

• A symmetric,  high-order  interpolation  is  used  obtain  spatial  distribution  of 
the  fine-block  variables  at  the  interface  in  order  to  prevent  artificial 
spurious  oscillations  in  the  flow  from  being  originated  from  the  interface. 

Studying  cases  show  that  the  multi-block  method  can  greatly  increase  the 
computational  efficiency  of  LBE  method  without  lost  of  accuracy. 

2.  Multi-time-relaxation  model  in  the  LBE  method 

In  computing  high  Re  flows,  the  computational  stability  will  become  a big 
issue.  Computational  stability  depends  on  the  LBE  models,  the  boundary 
conditions,  and  the  Re  number  of  flow.  In  the  present  study,  two  different 
LBE  models,  SRT  model  and  MRT  model,  are  studied  in  various  cases 
involving  singularity  at  high  Reynolds  number.  It  is  found  that 

• The  difference  between  the  two  models  can  be  non-local.  The  MRT 
model  in  general  provides  smoother  variations  of  the  macroscopic 
quantities  and  has  much  smaller  regions  of  the  oscillation  near  a 
singularity. 
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• Since  the  spatial  oscillation  is  often  accompanied  by  the  high  frequency 
pressure  (acoustic)  waves  in  transient  simulations,  the  MRT  model  also 
offers  a better  convergence  toward  steady  state  as  well. 

3.  Force  evaluation  on  the  solid  body  in  the  LBE  method 

Evaluating  the  force  acted  on  a body  by  the  fluid  is  a very  important  issue  in 
the  simulation.  In  conjunction  with  the  LBE  method,  two  methods  for  force 
evaluation  are  examined:  one  is  based  on  stress  integration  on  the  surface  and 
the  other  is  based  on  momentum  exchange  between  the  fluid  and  the  solid 
boundary.  The  momentum  exchange  method  is  found  to  be  very  simple  to 
implement  while  the  integration  of  stress  requires  tedious  evaluation  of  the 
details  of  the  surface  geometry  in  addition  to  the  extrapolation  for  stress 
related  variable  on  the  surface.  The  momentum  exchange  method  gives  a good 
result  and  thus  is  recommended  for  force  evaluation. 

4.  The  open  boundary  treatment 

An  open  boundary  condition  based  on  the  known  equilibrium  and  non- 
equilibrium values  was  proposed.  Compared  with  conventional  bounce 
scheme,  it  was  found  that  new  treatment  can 

• Avoid  strong  coupling  or  interaction  between  flow  variables  at  the 
boundary  and  at  the  interior  of  the  flow  field. 

• Improve  the  convergence  rate  to  steady  state,  stability  of  computation,  and 
quality  of  the  solution. 

5.  A unified  treatment  for  solid  wall  boundary 
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A unified  second-order  accurate  boundary  condition  treatment  for  the  solid 
wall  is  developed.  This  new  treatment  is: 

• Unified  for  A>0.5  and  A<0.5  and  there  is  no  discontinuity  in  the  boundary 
treatment. 

• Very  robust  to  implement  and  surfers  little  computational  instabilities 

• Overall  better  in  the  computational  performance  when  using  three-points 
interpolation  compared  with  other  boundary  treatments. 

6.  An  improved  moving  wall  boundary  condition 

Moving  boundary  problems  generate  additional  difficulties  during  the 
computation.,  which  include 

• Phase  change,  which  means  a node  in  the  fluid  may  become  a node  inside 
solid  and  vice  verse. 

• Change  of  the  Position  of  Singularity 

• The  change  of  the  number  of  nodes  which  are  used  in  the  force  evaluation 
using  momentum  exchange  method. 

A new  boundary  condition  specifically  designed  for  cases  with  moving  wall  is 
presented  to  solve  the  phase  change  problem.  As  the  boundary  crosses  the 
lattice  nodes,  analysis  shows  that  new  boundary  treatment  can 

• ensure  the  continuity  of  the  flow  variable  near  the  solid  boundary 

• ensure  the  continuity  of  momentum  exchange 

Computational  results  demonstrate  that  the  proposed  boundary  condition  can 
reduce  the  oscillation  caused  by  the  moving  boundary  when  compared  with 
the  solution  obtained  using  other  boundary  treatments. 
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For  the  moving  boundary  problem,  a stronger  interaction  exists  between  the  solid  and 
fluid  and  the  accuracy  of  the  force  acting  on  the  body  plays  even  a bigger  role  since  it 
feeds  back  to  the  flow  field  through  the  boundary  condition.  In  the  future  work,  more 
accurate  force  evaluation  scheme  needs  to  be  developed  to  account  for  the  changing 
configurations  of  the  solid  boundary.  Strategies  must  be  developed  to  reduce  the 
transient  oscillation  in  the  flow  variables  as  the  boundary  moves  in  the  flow  field. 
Especially,  there  is  noise  in  pressure  field  in  the  flow  near  a solid  body.  This  noise  exists 
even  in  the  flow  over  a stationary  body,  although  in  this  case  it  will  not  cause  big 
problem  in  simulation.  When  boundary  moves  this  kind  of  noise  will  generate  spatial  and 
temporal  pressure  oscillation  in  computation.  To  remove  such  kind  of  noise  in  the 
pressure  field  is  another  goal  in  the  study  of  moving  boundary  problem. 

In  the  present  study,  only  laminar  flows  are  being  considered.  When  Reynolds 
number  increase,  flows  will  become  turbulent,  so  it  is  necessary  to  introduce  turbulence 
model  in  computation  to  extend  the  application  of  the  LBE  method 
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